rm(list=ls())
gc()
setwd("/hpc/group/pbenfeylab/CheWei/CW_data/genesys/")

## Need seu4
suppressMessages(library(Seurat))
suppressMessages(library(cowplot))
suppressMessages(library(scattermore))
suppressMessages(library(scater))
suppressMessages(library(cowplot))
suppressMessages(library(RColorBrewer))
suppressMessages(library(grid))
suppressMessages(library(gplots))
suppressMessages(library(circular))
suppressMessages(library(ggplot2))
suppressMessages(library(ggnewscale))
suppressMessages(library(tidyverse))
suppressMessages(library(ComplexHeatmap))
suppressMessages(library(circlize))
suppressMessages(library(patchwork))

sessionInfo()

R version 4.2.2 (2022-10-31)
Platform: x86_64-conda-linux-gnu (64-bit)
Running under: AlmaLinux 9.3 (Shamrock Pampas Cat)

Matrix products: default
BLAS/LAPACK: /hpc/group/pbenfeylab/ch416/miniconda3/envs/seu4/lib/libopenblasp-r0.3.21.so

locale:
 [1] LC_CTYPE=en_US.UTF-8       LC_NUMERIC=C              
 [3] LC_TIME=en_US.UTF-8        LC_COLLATE=en_US.UTF-8    
 [5] LC_MONETARY=en_US.UTF-8    LC_MESSAGES=en_US.UTF-8   
 [7] LC_PAPER=en_US.UTF-8       LC_NAME=C                 
 [9] LC_ADDRESS=C               LC_TELEPHONE=C            
[11] LC_MEASUREMENT=en_US.UTF-8 LC_IDENTIFICATION=C       

attached base packages:
[1] grid      stats4    stats     graphics  grDevices utils     datasets 
[8] methods   base     

other attached packages:
 [1] patchwork_1.1.3             circlize_0.4.15            
 [3] ComplexHeatmap_2.14.0       forcats_0.5.2              
 [5] stringr_1.5.1               dplyr_1.1.3                
 [7] purrr_1.0.2                 readr_2.1.3                
 [9] tidyr_1.3.0                 tibble_3.2.1               
[11] tidyverse_1.3.2             ggnewscale_0.4.8           
[13] circular_0.4-95             gplots_3.1.3               
[15] RColorBrewer_1.1-3          scater_1.26.1              
[17] ggplot2_3.4.4               scuttle_1.8.0              
[19] SingleCellExperiment_1.20.0 SummarizedExperiment_1.28.0
[21] Biobase_2.58.0              GenomicRanges_1.50.0       
[23] GenomeInfoDb_1.34.8         IRanges_2.32.0             
[25] S4Vectors_0.36.0            BiocGenerics_0.44.0        
[27] MatrixGenerics_1.10.0       matrixStats_1.1.0          
[29] scattermore_1.2             cowplot_1.1.1              
[31] SeuratObject_4.1.3          Seurat_4.1.1.9001          

loaded via a namespace (and not attached):
  [1] utf8_1.2.4                spatstat.explore_3.2-5   
  [3] reticulate_1.34.0         tidyselect_1.2.0         
  [5] htmlwidgets_1.6.2         BiocParallel_1.32.5      
  [7] Rtsne_0.16                munsell_0.5.0            
  [9] ScaledMatrix_1.6.0        codetools_0.2-19         
 [11] ica_1.0-3                 pbdZMQ_0.3-8             
 [13] future_1.33.0             miniUI_0.1.1.1           
 [15] withr_2.5.2               spatstat.random_3.2-1    
 [17] colorspace_2.1-0          progressr_0.14.0         
 [19] uuid_1.1-0                ROCR_1.0-11              
 [21] tensor_1.5                listenv_0.9.0            
 [23] repr_1.1.4                GenomeInfoDbData_1.2.9   
 [25] polyclip_1.10-6           parallelly_1.36.0        
 [27] vctrs_0.6.4               generics_0.1.3           
 [29] timechange_0.1.1          doParallel_1.0.17        
 [31] R6_2.5.1                  clue_0.3-64              
 [33] ggbeeswarm_0.7.1          rsvd_1.0.5               
 [35] bitops_1.0-7              spatstat.utils_3.0-4     
 [37] DelayedArray_0.24.0       assertthat_0.2.1         
 [39] promises_1.2.1            scales_1.2.1             
 [41] googlesheets4_1.0.1       beeswarm_0.4.0           
 [43] gtable_0.3.4              beachmat_2.14.0          
 [45] globals_0.16.2            goftest_1.2-3            
 [47] rlang_1.1.2               GlobalOptions_0.1.2      
 [49] splines_4.2.2             lazyeval_0.2.2           
 [51] gargle_1.2.1              spatstat.geom_3.2-7      
 [53] broom_1.0.2               modelr_0.1.10            
 [55] reshape2_1.4.4            abind_1.4-5              
 [57] backports_1.4.1           httpuv_1.6.12            
 [59] tools_4.2.2               ellipsis_0.3.2           
 [61] ggridges_0.5.4            Rcpp_1.0.11              
 [63] plyr_1.8.9                base64enc_0.1-3          
 [65] sparseMatrixStats_1.10.0  zlibbioc_1.44.0          
 [67] RCurl_1.98-1.6            deldir_1.0-9             
 [69] GetoptLong_1.0.5          pbapply_1.7-2            
 [71] viridis_0.6.4             zoo_1.8-12               
 [73] haven_2.5.1               ggrepel_0.9.4            
 [75] cluster_2.1.4             fs_1.6.3                 
 [77] magrittr_2.0.3            data.table_1.14.8        
 [79] RSpectra_0.16-1           reprex_2.0.2             
 [81] lmtest_0.9-40             RANN_2.6.1               
 [83] googledrive_2.0.0         mvtnorm_1.1-3            
 [85] fitdistrplus_1.1-11       hms_1.1.2                
 [87] mime_0.12                 evaluate_0.23            
 [89] xtable_1.8-4              readxl_1.4.1             
 [91] shape_1.4.6               fastDummies_1.7.3        
 [93] gridExtra_2.3             compiler_4.2.2           
 [95] KernSmooth_2.23-20        crayon_1.5.2             
 [97] htmltools_0.5.7           tzdb_0.3.0               
 [99] later_1.3.1               lubridate_1.9.0          
[101] DBI_1.1.3                 dbplyr_2.2.1             
[103] MASS_7.3-58.3             boot_1.3-28.1            
[105] Matrix_1.5-4              cli_3.6.1                
[107] parallel_4.2.2            igraph_1.5.1             
[109] pkgconfig_2.0.3           sp_2.1-1                 
[111] IRdisplay_1.1             plotly_4.10.3            
[113] spatstat.sparse_3.0-3     foreach_1.5.2            
[115] xml2_1.3.3                vipor_0.4.5              
[117] XVector_0.38.0            rvest_1.0.3              
[119] digest_0.6.33             sctransform_0.4.1        
[121] RcppAnnoy_0.0.21          spatstat.data_3.0-3      
[123] cellranger_1.1.0          leiden_0.4.3             
[125] uwot_0.1.16               DelayedMatrixStats_1.20.0
[127] shiny_1.7.5.1             gtools_3.9.4             
[129] rjson_0.2.21              lifecycle_1.0.4          
[131] nlme_3.1-162              jsonlite_1.8.7           
[133] BiocNeighbors_1.16.0      viridisLite_0.4.2        
[135] fansi_1.0.5               pillar_1.9.0             
[137] lattice_0.21-8            fastmap_1.1.1            
[139] httr_1.4.7                survival_3.4-0           
[141] glue_1.6.2                iterators_1.0.14         
[143] png_0.1-8                 stringi_1.8.1            
[145] RcppHNSW_0.5.0            BiocSingular_1.14.0      
[147] caTools_1.18.2            IRkernel_1.3.1.9000      
[149] irlba_2.3.5.1             future.apply_1.11.0

wanted_TFs <- read.csv("./Kay_TF_thalemine_annotations.csv")

nrow(wanted_TFs)

## Make TF names unique
wanted_TFs$Name[wanted_TFs$GeneID=="AT2G33880"]="WOX9"
wanted_TFs$Name[wanted_TFs$GeneID=="AT2G45160"]="SCL27"
wanted_TFs$Name[wanted_TFs$GeneID=="AT5G04410"]="NAC78"
wanted_TFs$Name[wanted_TFs$GeneID=="AT3G29035"]="ORS1"
wanted_TFs$Name[wanted_TFs$GeneID=="AT2G02540"]="ZHD3"
wanted_TFs$Name[wanted_TFs$GeneID=="AT3G16500"]="IAA26"
wanted_TFs$Name[wanted_TFs$GeneID=="AT5G09740"]="HAG5"
wanted_TFs$Name[wanted_TFs$GeneID=="AT4G24660"]="ZHD2"
wanted_TFs$Name[wanted_TFs$GeneID=="AT5G46880"]="HDG5"
wanted_TFs$Name[wanted_TFs$GeneID=="AT1G28420"]="RLT1"
wanted_TFs$Name[wanted_TFs$GeneID=="AT1G14580"]="BLJ"
wanted_TFs$Name[wanted_TFs$GeneID=="AT3G45260"]="BIB"
wanted_TFs$Name[wanted_TFs$GeneID=="AT2G02070"]="RVN"
wanted_TFs$Name[wanted_TFs$GeneID=="AT2G28160"]="FIT"
wanted_TFs$Name[wanted_TFs$GeneID=="AT1G68360"]="GIS3"
wanted_TFs$Name[wanted_TFs$GeneID=="AT1G20640"]="NLP4"
wanted_TFs$Name[wanted_TFs$GeneID=="AT5G05550"]="VFP5"
wanted_TFs$Name[wanted_TFs$GeneID=="AT3G59470"]="FRF1"
wanted_TFs$Name[wanted_TFs$GeneID=="AT5G15150"]="HAT7"
wanted_TFs$Name[wanted_TFs$GeneID=="AT5G14750"]="WER"
wanted_TFs$Name[wanted_TFs$GeneID=="AT1G75710"]="BRON"
wanted_TFs$Name[wanted_TFs$GeneID=="AT1G74500"]="TMO7"
wanted_TFs$Name[wanted_TFs$GeneID=="AT2G12646"]="RITF1"
wanted_TFs$Name[wanted_TFs$GeneID=="AT3G48100"]="ARR5"
wanted_TFs$Name[wanted_TFs$GeneID=="AT4G16141"]="GATA17L"
wanted_TFs$Name[wanted_TFs$GeneID=="AT5G65640"]="NFL"
wanted_TFs$Name[wanted_TFs$GeneID=="AT1G62700"]="VND5"
wanted_TFs$Name[wanted_TFs$GeneID=="AT4G36160"]="VND2"
wanted_TFs$Name[wanted_TFs$GeneID=="AT5G66300"]="VND3"
wanted_TFs$Name[wanted_TFs$GeneID=="AT1G12260"]="VND4"
wanted_TFs$Name[wanted_TFs$GeneID=="AT5G62380"]="VND6"

## TTG1
wanted_TFs$Name[wanted_TFs$GeneID=="AT5G24520"]

## SCRAMBLED
wanted_TFs$Name[wanted_TFs$GeneID=="AT1G11130"]

## CAPRICE
wanted_TFs$Name[wanted_TFs$GeneID=="AT2G46410"]

stem2pro <- read.csv("./TF_GRN_centrality_t0-t1_zscore3.csv")
pro2trans <- read.csv("./TF_GRN_centrality_t1-t3_zscore3.csv")
trans2el <- read.csv("./TF_GRN_centrality_t3-t5_zscore3.csv")
el2el <- read.csv("./TF_GRN_centrality_t5-t7_zscore3.csv")
el2mat <- read.csv("./TF_GRN_centrality_t7-t9_zscore3.csv")

head(stem2pro)

head(pro2trans)

head(trans2el)

head(el2el)

head(el2mat)

head(trans2el)

min_max_normalize <- function(data) {
  min_val <- min(data)
  max_val <- max(data)
  normalized_data <- (data - min_val) / (max_val - min_val)
  return(normalized_data)
}

summary(stem2pro$tri_betweenness_centrality)

   Min. 1st Qu.  Median    Mean 3rd Qu.    Max. 
0.00000 0.00000 0.00000 0.03773 0.00000 0.98370

summary(min_max_normalize(stem2pro$tri_betweenness_centrality))

   Min. 1st Qu.  Median    Mean 3rd Qu.    Max. 
0.00000 0.00000 0.00000 0.03835 0.00000 1.00000

ncol(stem2pro)

stem2pro <- as.data.frame(cbind(stem2pro[,1],apply(stem2pro[,3:62],2,min_max_normalize)))
pro2trans <- as.data.frame(cbind(pro2trans[,1],apply(pro2trans[,3:62],2,min_max_normalize)))
trans2el <- as.data.frame(cbind(trans2el[,1],apply(trans2el[,3:62],2,min_max_normalize)))
el2el <- as.data.frame(cbind(el2el[,1],apply(el2el[,3:62],2,min_max_normalize)))
el2mat <- as.data.frame(cbind(el2mat[,1],apply(el2mat[,3:62],2,min_max_normalize)))

dat <- stem2pro %>%
full_join(pro2trans, by = "V1") %>%
full_join(trans2el, by = "V1") %>%
full_join(el2el, by = "V1") %>%
full_join(el2mat, by = "V1")

dat[is.na(dat)] <- 0

n <- c('atri_degree_centrality','atri_out_centrality','atri_in_centrality','atri_betweenness_centrality','atri_closeness_centrality','atri_eigenvector_centrality',
       'tri_degree_centrality','tri_out_centrality','tri_in_centrality','tri_betweenness_centrality','tri_closeness_centrality','tri_eigenvector_centrality',
       'lrc_degree_centrality','lrc_out_centrality','lrc_in_centrality','lrc_betweenness_centrality','lrc_closeness_centrality','lrc_eigenvector_centrality',
       'cor_degree_centrality','cor_out_centrality','cor_in_centrality','cor_betweenness_centrality','cor_closeness_centrality','cor_eigenvector_centrality',
       'end_degree_centrality','end_out_centrality','end_in_centrality','end_betweenness_centrality','end_closeness_centrality','end_eigenvector_centrality',
       'per_degree_centrality','per_out_centrality','per_in_centrality','per_betweenness_centrality','per_closeness_centrality','per_eigenvector_centrality',
       'pro_degree_centrality','pro_out_centrality','pro_in_centrality','pro_betweenness_centrality','pro_closeness_centrality','pro_eigenvector_centrality',
       'xyl_degree_centrality','xyl_out_centrality','xyl_in_centrality','xyl_betweenness_centrality','xyl_closeness_centrality','xyl_eigenvector_centrality',
       'phl_degree_centrality','phl_out_centrality','phl_in_centrality','phl_betweenness_centrality','phl_closeness_centrality','phl_eigenvector_centrality',
       'col_degree_centrality','col_out_centrality','col_in_centrality','col_betweenness_centrality','col_closeness_centrality','col_eigenvector_centrality')

colnames(dat) <- c("TF",gsub("$","_1",n), gsub("$","_2",n),gsub("$","_3",n),gsub("$","_4",n),gsub("$","_5",n))

GeneID <- wanted_TFs$GeneID[match(dat$TF, wanted_TFs$Name)]

dat <- cbind(GeneID, dat)

head(dat)

numz <- function(x){
    sum(x==0)/length(x)
}

dat$combined_score <- min_max_normalize(rowSums(apply(dat[,grep("centrality",colnames(dat))],2,as.numeric)))
dat$celltype_specificity <- min_max_normalize(apply(apply(dat[,grep("centrality",colnames(dat))],2,as.numeric),1,numz))
dat$weighted_score <- dat$combined_score + dat$celltype_specificity
dat <- dat %>% arrange(desc(weighted_score))

head(dat)

write.csv(dat,"TF_GRN_centrality_all_transition_zscore3.csv", quote=FALSE, row.names=FALSE)

dat_genesys <- dat

stem2pro <- read.csv("./scRNAseq_TF_GRN_centrality_t0-t1_zscore3.csv")
pro2trans <- read.csv("./scRNAseq_TF_GRN_centrality_t1-t3_zscore3.csv")
trans2el <- read.csv("./scRNAseq_TF_GRN_centrality_t3-t5_zscore3.csv")
el2el <- read.csv("./scRNAseq_TF_GRN_centrality_t5-t7_zscore3.csv")
el2mat <- read.csv("./scRNAseq_TF_GRN_centrality_t7-t9_zscore3.csv")

head(stem2pro)

head(pro2trans)

head(trans2el)

head(el2el)

head(el2mat)

stem2pro <- as.data.frame(cbind(stem2pro[,1],apply(stem2pro[,3:62],2,min_max_normalize)))
pro2trans <- as.data.frame(cbind(pro2trans[,1],apply(pro2trans[,3:62],2,min_max_normalize)))
trans2el <- as.data.frame(cbind(trans2el[,1],apply(trans2el[,3:62],2,min_max_normalize)))
el2el <- as.data.frame(cbind(el2el[,1],apply(el2el[,3:62],2,min_max_normalize)))
el2mat <- as.data.frame(cbind(el2mat[,1],apply(el2mat[,3:62],2,min_max_normalize)))

dat <- stem2pro %>%
full_join(pro2trans, by = "V1") %>%
full_join(trans2el, by = "V1") %>%
full_join(el2el, by = "V1") %>%
full_join(el2mat, by = "V1")

dat[is.na(dat)] <- 0

n <- c('atri_degree_centrality','atri_out_centrality','atri_in_centrality','atri_betweenness_centrality','atri_closeness_centrality','atri_eigenvector_centrality',
       'tri_degree_centrality','tri_out_centrality','tri_in_centrality','tri_betweenness_centrality','tri_closeness_centrality','tri_eigenvector_centrality',
       'lrc_degree_centrality','lrc_out_centrality','lrc_in_centrality','lrc_betweenness_centrality','lrc_closeness_centrality','lrc_eigenvector_centrality',
       'cor_degree_centrality','cor_out_centrality','cor_in_centrality','cor_betweenness_centrality','cor_closeness_centrality','cor_eigenvector_centrality',
       'end_degree_centrality','end_out_centrality','end_in_centrality','end_betweenness_centrality','end_closeness_centrality','end_eigenvector_centrality',
       'per_degree_centrality','per_out_centrality','per_in_centrality','per_betweenness_centrality','per_closeness_centrality','per_eigenvector_centrality',
       'pro_degree_centrality','pro_out_centrality','pro_in_centrality','pro_betweenness_centrality','pro_closeness_centrality','pro_eigenvector_centrality',
       'xyl_degree_centrality','xyl_out_centrality','xyl_in_centrality','xyl_betweenness_centrality','xyl_closeness_centrality','xyl_eigenvector_centrality',
       'phl_degree_centrality','phl_out_centrality','phl_in_centrality','phl_betweenness_centrality','phl_closeness_centrality','phl_eigenvector_centrality',
       'col_degree_centrality','col_out_centrality','col_in_centrality','col_betweenness_centrality','col_closeness_centrality','col_eigenvector_centrality')

colnames(dat) <- c("TF",gsub("$","_1",n), gsub("$","_2",n),gsub("$","_3",n),gsub("$","_4",n),gsub("$","_5",n))

GeneID <- wanted_TFs$GeneID[match(dat$TF, wanted_TFs$Name)]

dat <- cbind(GeneID, dat)

head(dat)

dat$combined_score <- min_max_normalize(rowSums(apply(dat[,grep("centrality",colnames(dat))],2,as.numeric)))
dat$celltype_specificity <- min_max_normalize(apply(apply(dat[,grep("centrality",colnames(dat))],2,as.numeric),1,numz))
dat$weighted_score <- dat$combined_score + dat$celltype_specificity
dat <- dat %>% arrange(desc(weighted_score))

head(dat)

write.csv(dat,"scRNAseq_TF_GRN_centrality_all_transition_zscore3.csv", quote=FALSE, row.names=FALSE)

dat_raw <- dat

gene_list <- read.table('./gene_list_1108.csv', sep=",", header = TRUE)

exptf <- intersect(gene_list$features, wanted_TFs$GeneID)
length(exptf)

sfun <- read.csv('string_functional_annotations.tsv', sep="\t", header=TRUE)
sann <- read.csv('string_protein_annotations.tsv', sep="\t", header=TRUE)

gsgo1 <- unique(sfun[grep('root|xylem|phloem|procambium|pericycle|vascular|vasculature|stele|tracheary|sieve|trichoblast|atrichoblast|epidermis|epidermal tissue|lateral root cap|root hair|trichome|cortex|endodermis|ground tissue|columella|quiescent center',sfun$term.description, ignore.case=TRUE),]$X.node)
gsgo2 <- unique(sann[grep('root|xylem|phloem|procambium|pericycle|vascular|vasculature|stele|tracheary|sieve|trichoblast|atrichoblast|epidermis|epidermal tissue|lateral root cap|root hair|trichome|cortex|endodermis|ground tissue|columella|quiescent center',sann$domain_summary_url, ignore.case=TRUE),]$X.node)
gsgo <- sort(unique(c(gsgo1, gsgo2)))

#write.csv(data.frame(GeneID=gsgo),"./Gold_Standard_Root_TF_StringDB.csv", quote=FALSE, row.names=FALSE)

length(gsgo)
gsgo

gsgo <- gsub(",.*$","",gsub("^.*ath:","",sann[match(gsgo, sann$X.node),]$other_names_and_aliases))
gsgo[which(gsgo=='831248')]='AT5G14000'
gsgo

sann[which(sann$X.node=="HAT7"),]

sann[grep("831248",sann$other_names_and_aliases, ignore.case = TRUE),]

sfun <- read.csv('./functional_GS_root_unique.csv', sep=",", header=TRUE)

gsgo_miniex <- intersect(exptf, sfun$GeneID)
length(gsgo_miniex)

gsgo <- intersect(gsgo, gsgo_miniex)

length(gsgo)

gsgo

write.csv(data.frame(GeneID=gsgo),"./Gold_Standard_Root_TF_MINI_EX_StringDB.csv", quote=FALSE, row.names=FALSE)

r50 <- 70
numz <- function(x){
    sum(x==0)/length(x)
}

genesys <- dat_genesys

run_r50_genesys <- function(x){
    genesys$combined_score <- min_max_normalize(rowSums(apply(genesys[,grep(x,colnames(genesys))],2,as.numeric)))
    genesys <- genesys %>% arrange(desc(combined_score))
    count <- 0
    for (i in seq(nrow(genesys))){
        if (genesys$GeneID[i] %in% gsgo){
            count <- count +1
            if (count == r50){
                return(i)
                break
            }
        }
    }
}

## all cell type : in centrality + celltype & dev stage specificity
#dat$combined_score <- min_max_normalize(rowSums(apply(dat[,grep("in_centrality",colnames(dat))],2,as.numeric)))
#dat$celltype_specificity <- min_max_normalize(apply(apply(dat[,grep("in_centrality",colnames(dat))],2,as.numeric),1,numz))
#dat$weighted_score <- dat$combined_score + dat$celltype_specificity
#dat <- dat %>% arrange(desc(weighted_score))

#count <- 0
#for (i in seq(nrow(dat))){
#    if (dat$GeneID[i] %in% gsgo){
#        count <- count +1
#        if (count == r50){
#            print(i)
#            break
#        }
#    }
#}

raw <- dat_raw

run_r50_scRNAseq <- function(x){
    raw$combined_score <- min_max_normalize(rowSums(apply(raw[,grep(x,colnames(raw))],2,as.numeric)))
    raw <- raw %>% arrange(desc(combined_score))
    count <- 0
    for (i in seq(nrow(raw))){
        if (raw$GeneID[i] %in% gsgo){
            count <- count +1
            if (count == r50){
                return(i)
                break
            }
        }
    }
}

prepros <- function(x){
    dat <- read.csv(x)
    #dat <- dat %>% filter(role=="Connector Hub") 
    dat <- dat[,c(1, grep("centrality",colnames(dat)), 11)]
    return(dat)
}

atri <- prepros("../celloracle/atrichoblast_X_root_integrated_celloracle_gene_score_iGRN.csv")
tri <- prepros("../celloracle/trichoblast_X_root_integrated_celloracle_gene_score_iGRN.csv")
lrc <- prepros("../celloracle/lrc_X_root_integrated_celloracle_gene_score_iGRN.csv")
cor <- prepros("../celloracle/cortex_X_root_integrated_celloracle_gene_score_iGRN.csv")
end <- prepros("../celloracle/endodermis_X_root_integrated_celloracle_gene_score_iGRN.csv")
per <- prepros("../celloracle/pericycle_X_root_integrated_celloracle_gene_score_iGRN.csv")
pro <- prepros("../celloracle/procambium_X_root_integrated_celloracle_gene_score_iGRN.csv")
xyl <- prepros("../celloracle/xylem_X_root_integrated_celloracle_gene_score_iGRN.csv")
phl <- prepros("../celloracle/phloem_X_root_integrated_celloracle_gene_score_iGRN.csv")
col <- prepros("../celloracle/columella_X_root_integrated_celloracle_gene_score_iGRN.csv")

dat <- rbind(atri, tri, lrc, cor, end, per, pro, xyl, phl, col)
## Remove those not present in the dataset
dat <- dat[!is.na(match(dat$X,exptf)),]
dat <- dat %>% group_by(cluster) %>% reframe(GeneID=X,degree_centrality = min_max_normalize(degree_centrality_all), 
                                  in_centrality = min_max_normalize(degree_centrality_in), 
                                  out_centrality = min_max_normalize(degree_centrality_out),
                                 betweenness_centrality = min_max_normalize(betweenness_centrality),
                                 #closeness_centrality = min_max_normalize(closeness_centrality),
                                 eigenvector_centrality = min_max_normalize(eigenvector_centrality))
dat <- dat %>% group_by(GeneID) %>% reframe(degree_centrality = sum(degree_centrality), 
                                  in_centrality = sum(in_centrality), 
                                  out_centrality = sum(out_centrality),
                                 betweenness_centrality = sum(betweenness_centrality),
                                 #closeness_centrality = sum(closeness_centrality),
                                 eigenvector_centrality = sum(eigenvector_centrality))
dat <- cbind(dat[,1],apply(dat[,2:6],2,min_max_normalize))

head(dat)

celloracle_scRNAseq <- dat

atri <- prepros("../celloracle/atrichoblast_P_root_integrated_celloracle_gene_score_iGRN.csv")
tri <- prepros("../celloracle/trichoblast_P_root_integrated_celloracle_gene_score_iGRN.csv")
lrc <- prepros("../celloracle/lrc_P_root_integrated_celloracle_gene_score_iGRN.csv")
cor <- prepros("../celloracle/cortex_P_root_integrated_celloracle_gene_score_iGRN.csv")
end <- prepros("../celloracle/endodermis_P_root_integrated_celloracle_gene_score_iGRN.csv")
per <- prepros("../celloracle/pericycle_P_root_integrated_celloracle_gene_score_iGRN.csv")
pro <- prepros("../celloracle/procambium_P_root_integrated_celloracle_gene_score_iGRN.csv")
xyl <- prepros("../celloracle/xylem_P_root_integrated_celloracle_gene_score_iGRN.csv")
phl <- prepros("../celloracle/phloem_P_root_integrated_celloracle_gene_score_iGRN.csv")
col <- prepros("../celloracle/columella_P_root_integrated_celloracle_gene_score_iGRN.csv")

dat <- rbind(atri, tri, lrc, cor, end, per, pro, xyl, phl, col)
## Remove those not present in the dataset
dat <- dat[!is.na(match(dat$X,exptf)),]
dat <- dat %>% group_by(cluster) %>% reframe(GeneID=X,degree_centrality = min_max_normalize(degree_centrality_all), 
                                  in_centrality = min_max_normalize(degree_centrality_in), 
                                  out_centrality = min_max_normalize(degree_centrality_out),
                                 betweenness_centrality = min_max_normalize(betweenness_centrality),
                                 #closeness_centrality = min_max_normalize(closeness_centrality),
                                 eigenvector_centrality = min_max_normalize(eigenvector_centrality))
dat <- dat %>% group_by(GeneID) %>% reframe(degree_centrality = sum(degree_centrality), 
                                  in_centrality = sum(in_centrality), 
                                  out_centrality = sum(out_centrality),
                                 betweenness_centrality = sum(betweenness_centrality),
                                 #closeness_centrality = sum(closeness_centrality),
                                 eigenvector_centrality = sum(eigenvector_centrality))
dat <- cbind(dat[,1],apply(dat[,2:6],2,min_max_normalize))

head(dat)

celloracle_GeneSys <- dat

run_r50_celloracle <- function(x, dat){
    celloracle <- dat
    celloracle$combined_score <- celloracle[,grep(x,colnames(celloracle))]
    celloracle <- celloracle %>% arrange(desc(combined_score))
    count <- 0
    for (i in seq(nrow(celloracle))){
        if (celloracle$GeneID[i] %in% gsgo){
            count <- count +1
            if (count == r50){
                return(i)
                break
            }
        }
    }
}

de_scRNAseq <- read.csv("Input_X_root_markers_ROC.csv")

de_scRNAseq <- de_scRNAseq %>% arrange(pct.diff_rank) %>% arrange(avg_diff_rank)%>% arrange(myAUC_rank)%>% arrange(combined_rank)

## Remove those not present in the dataset
de_scRNAseq <- de_scRNAseq[!is.na(match(de_scRNAseq$gene,exptf)),]

dat <- de_scRNAseq %>% group_by(gene) %>% reframe(myAUC_rank = mean(myAUC_rank),avg_diff_rank = mean(avg_diff_rank))
dat <- dat %>% arrange(myAUC_rank)
colnames(dat) <- c("GeneID","myAUC_rank","avg_diff_rank")

de_scRNAseq <- dat

de <- read.csv("GeneSys_P_root_markers_ROC.csv")

de <- de %>% arrange(pct.diff_rank) %>% arrange(avg_diff_rank)%>% arrange(myAUC_rank)%>% arrange(combined_rank)

## Remove those not present in the dataset
de <- de[!is.na(match(de$gene,exptf)),]

head(de)

dat <- de %>% group_by(gene) %>% reframe(myAUC_rank = mean(myAUC_rank),avg_diff_rank = mean(avg_diff_rank))
dat <- dat %>% arrange(myAUC_rank)
colnames(dat) <- c("GeneID","myAUC_rank","avg_diff_rank")

de_GeneSys <- dat

run_r50_de <- function(x, dat){
    de <- dat
    de$combined_score <- de[,grep(x,colnames(de))]
    de <- de %>% arrange(desc(combined_score))
    count <- 0
    for (i in seq(nrow(de))){
        if (de$GeneID[i] %in% gsgo){
            count <- count +1
            if (count == r50){
                return(i)
                break
            }
        }
    }
}

## Expressed TFs permutation
R50_permutation <- c()
for (j in 1:1000){
ran <- sample(exptf, length(exptf), replace=FALSE)
count <- 0
for (i in seq(length(ran))){
    if (ran[i] %in% gsgo){
        count <- count +1
        if (count == r50){
            R50_permutation <- c(R50_permutation,i)
            break
        }
    }
}
}

toplt <- data.frame(Methods=c("LIMA degree centrality (GeneSys)", "LIMA out centrality (GeneSys)", "LIMA in centrality (GeneSys)", "LIMA betweenness centrality (GeneSys)", 
                 "LIMA eigenvector centrality (GeneSys)","LIMA degree centrality (scRNAseq)", "LIMA out centrality (scRNAseq)", "LIMA in centrality (scRNAseq)", "LIMA betweenness centrality (scRNAseq)", 
                "LIMA eigenvector centrality (scRNAseq)", "CellOracle degree centrality (scRNAseq)", "CellOracle out centrality (scRNAseq)", 
                 "CellOracle in centrality (scRNAseq)", "CellOracle betweenness centrality (scRNAseq)", "CellOracle eigenvector centrality (scRNAseq)",
                 "CellOracle degree centrality (GeneSys)", "CellOracle out centrality (GeneSys)", 
                 "CellOracle in centrality (GeneSys)", "CellOracle betweenness centrality (GeneSys)", "CellOracle eigenvector centrality (GeneSys)",
                 "DE myAUC rank (scRNAseq)", "DE avg diff rank (scRNAseq)", "DE myAUC rank (GeneSys)", "DE avg diff rank (GeneSys)", "Permutation"),
           R50=c(run_r50_genesys('degree_centrality'),run_r50_genesys('out_centrality'),run_r50_genesys('in_centrality'),
                 run_r50_genesys('betweenness_centrality'),run_r50_genesys('eigenvector_centrality'),
                 run_r50_scRNAseq('degree_centrality'),run_r50_scRNAseq('out_centrality'),run_r50_scRNAseq('in_centrality'),
                 run_r50_scRNAseq('betweenness_centrality'),run_r50_scRNAseq('eigenvector_centrality'),
                 run_r50_celloracle('degree_centrality', celloracle_scRNAseq),run_r50_celloracle('out_centrality', celloracle_scRNAseq),run_r50_celloracle('in_centrality', celloracle_scRNAseq),
                 run_r50_celloracle('betweenness_centrality', celloracle_scRNAseq),run_r50_celloracle('eigenvector_centrality', celloracle_scRNAseq),
                 run_r50_celloracle('degree_centrality', celloracle_GeneSys),run_r50_celloracle('out_centrality', celloracle_GeneSys),run_r50_celloracle('in_centrality', celloracle_GeneSys),
                 run_r50_celloracle('betweenness_centrality', celloracle_GeneSys),run_r50_celloracle('eigenvector_centrality', celloracle_GeneSys),
                 run_r50_de('myAUC_rank', de_scRNAseq),run_r50_de('avg_diff_rank', de_scRNAseq),run_r50_de('myAUC_rank', de_GeneSys),run_r50_de('avg_diff_rank', de_GeneSys), mean(R50_permutation)))

toplt

options(repr.plot.width=10, repr.plot.height=8)
ggplot(toplt, aes(x=reorder(Methods, R50, decreasing = TRUE),  y=R50)) + geom_point(size=4)+ 
  labs(title="TF Prioritization Performance",x="", y = "R50")+
  theme_classic()+
  theme(axis.text.x = element_text(color = "black", 
                           size = 16, angle = 0),
          axis.text.y = element_text( color = "black", 
                           size = 16, angle = 0),
        plot.title = element_text(color="black", size=18, face="bold"),
        axis.title.x = element_text(color="black", size=16, face="bold"))+ coord_flip()

toplt <- data.frame(Methods=factor(c(rep("LIMA (GeneSys)",5), rep("LIMA (scRNAseq)",5), rep("CellOracle (scRNAseq)",5), rep("CellOracle (GeneSys)",5), rep("DE (scRNAseq)",2), rep("DE (GeneSys)",2), "Permutation"), 
                                   levels=c("Permutation", "CellOracle (scRNAseq)", "LIMA (scRNAseq)", "CellOracle (GeneSys)", "DE (scRNAseq)", "LIMA (GeneSys)", "DE (GeneSys)")), R50=toplt$R50)

options(repr.plot.width=6, repr.plot.height=5)
ggplot(toplt, aes(x=Methods, y=R50))+
  geom_boxplot(width=0.4)+ 
  labs(title="TF Prioritization Performance",x="", y = "R50")+
  theme_classic()+
  theme(axis.text.x = element_text(color = "black", 
                           size = 16, angle = 0),
          axis.text.y = element_text( color = "black", 
                           size = 16, angle = 0),
        plot.title = element_text(color="black", size=18, face="bold"),
        axis.title.x = element_text(color="black", size=16, face="bold"))+ coord_flip()

toplt <- data.frame(Methods=factor(c(rep("LIMA (GeneSys)",5), rep("LIMA (scRNAseq)",5), rep("CellOracle (scRNAseq)",5), rep("CellOracle (GeneSys)",5), rep("DE (scRNAseq)",2), rep("DE (GeneSys)",2), "Permutation"), 
                                   levels=c("Permutation", "CellOracle (scRNAseq)", "CellOracle (GeneSys)", "LIMA (scRNAseq)", "LIMA (GeneSys)", "DE (scRNAseq)", "DE (GeneSys)")), R50=toplt$R50)

options(repr.plot.width=6, repr.plot.height=5)
ggplot(toplt, aes(x=Methods, y=R50))+
  geom_boxplot(width=0.4)+ 
  labs(title="TF Prioritization Performance",x="", y = "R50")+
  theme_classic()+
  theme(axis.text.x = element_text(color = "black", 
                           size = 16, angle = 0),
          axis.text.y = element_text( color = "black", 
                           size = 16, angle = 0),
        plot.title = element_text(color="black", size=18, face="bold"),
        axis.title.x = element_text(color="black", size=16, face="bold"))+ coord_flip()

toplt <- data.frame(Methods=factor(c(rep("LIMA (GeneSys)",5), rep("LIMA (scRNAseq)",5), rep("CellOracle (scRNAseq)",5), rep("CellOracle (GeneSys)",5), rep("DE (scRNAseq)",2), rep("DE (GeneSys)",2), "Permutation"), 
                                   levels=c("Permutation", "CellOracle (scRNAseq)", "LIMA (scRNAseq)", "DE (scRNAseq)", "CellOracle (GeneSys)", "LIMA (GeneSys)", "DE (GeneSys)")), R50=toplt$R50)

options(repr.plot.width=6, repr.plot.height=5)
ggplot(toplt, aes(x=Methods, y=R50))+
  geom_boxplot(width=0.4)+ 
  labs(title="TF Prioritization Performance",x="", y = "R50")+
  theme_classic()+
  theme(axis.text.x = element_text(color = "black", 
                           size = 16, angle = 0),
          axis.text.y = element_text( color = "black", 
                           size = 16, angle = 0),
        plot.title = element_text(color="black", size=18, face="bold"),
        axis.title.x = element_text(color="black", size=16, face="bold"))+ coord_flip()

gsgo_miniex_stringdb <- gsgo

gsgo <- read.csv("./Gold_Standard_Root_TF_StringDB.csv", header=TRUE)
gsgo <- gsgo[which(gsgo[,5]=="Yes"),]
gsgo$preferredName
gsgo <- gsgo$GeneID
gsgo <- intersect(gsgo, gsgo_miniex_stringdb)

length(gsgo)

r50 <- 16

## Define genesys
run_r50_genesys <- function(x){
    genesys$ct_score <- min_max_normalize(rowSums(apply(genesys[,grep(paste0('xyl_',x,'|phl_',x,'|pro_',x),colnames(genesys))],2,as.numeric)))
    genesys$combined_score <- genesys$ct_score
    
    genesys <- genesys %>% arrange(desc(combined_score))
    count <- 0
    for (i in seq(nrow(genesys))){
        if (genesys$GeneID[i] %in% gsgo){
            count <- count +1
            if (count == r50){
                return(i)
                break
            }
        }
    }
}

## Define raw
run_r50_scRNAseq <- function(x){
    raw$ct_score <- min_max_normalize(rowSums(apply(raw[,grep(paste0('xyl_',x,'|phl_',x,'|pro_',x),colnames(raw))],2,as.numeric)))
    raw$combined_score <- raw$ct_score
    
    raw <- raw %>% arrange(desc(combined_score))
    count <- 0
    for (i in seq(nrow(raw))){
        if (raw$GeneID[i] %in% gsgo){
            count <- count +1
            if (count == r50){
                return(i)
                break
            }
        }
    }
}

## Define celloracle
pro <- prepros("../celloracle/procambium_X_root_integrated_celloracle_gene_score_iGRN.csv")
xyl <- prepros("../celloracle/xylem_X_root_integrated_celloracle_gene_score_iGRN.csv")
phl <- prepros("../celloracle/phloem_X_root_integrated_celloracle_gene_score_iGRN.csv")
## Define celloracle
dat <- rbind(pro, xyl, phl)
## Remove those not present in the dataset
dat <- dat[!is.na(match(dat$X,exptf)),]
dat <- dat %>% group_by(cluster) %>% reframe(GeneID=X,degree_centrality = min_max_normalize(degree_centrality_all), 
                                  in_centrality = min_max_normalize(degree_centrality_in), 
                                  out_centrality = min_max_normalize(degree_centrality_out),
                                 betweenness_centrality = min_max_normalize(betweenness_centrality),
                                 eigenvector_centrality = min_max_normalize(eigenvector_centrality))
dat <- dat %>% group_by(GeneID) %>% reframe(degree_centrality = sum(degree_centrality), 
                                  in_centrality = sum(in_centrality), 
                                  out_centrality = sum(out_centrality),
                                 betweenness_centrality = sum(betweenness_centrality),
                                 eigenvector_centrality = sum(eigenvector_centrality))
dat <- cbind(dat[,1],apply(dat[,2:6],2,min_max_normalize))

celloracle_scRNAseq <- dat

## Define celloracle
pro <- prepros("../celloracle/procambium_P_root_integrated_celloracle_gene_score_iGRN.csv")
xyl <- prepros("../celloracle/xylem_P_root_integrated_celloracle_gene_score_iGRN.csv")
phl <- prepros("../celloracle/phloem_P_root_integrated_celloracle_gene_score_iGRN.csv")
## Define celloracle
dat <- rbind(pro, xyl, phl)
## Remove those not present in the dataset
dat <- dat[!is.na(match(dat$X,exptf)),]
dat <- dat %>% group_by(cluster) %>% reframe(GeneID=X,degree_centrality = min_max_normalize(degree_centrality_all), 
                                  in_centrality = min_max_normalize(degree_centrality_in), 
                                  out_centrality = min_max_normalize(degree_centrality_out),
                                 betweenness_centrality = min_max_normalize(betweenness_centrality),
                                 eigenvector_centrality = min_max_normalize(eigenvector_centrality))
dat <- dat %>% group_by(GeneID) %>% reframe(degree_centrality = sum(degree_centrality), 
                                  in_centrality = sum(in_centrality), 
                                  out_centrality = sum(out_centrality),
                                 betweenness_centrality = sum(betweenness_centrality),
                                 eigenvector_centrality = sum(eigenvector_centrality))
dat <- cbind(dat[,1],apply(dat[,2:6],2,min_max_normalize))

celloracle_GeneSys <- dat

run_r50_celloracle <- function(x, dat){
    celloracle <- dat
    celloracle$combined_score <- celloracle[,grep(x,colnames(celloracle))]
    celloracle <- celloracle %>% arrange(desc(combined_score))
    count <- 0
    for (i in seq(nrow(celloracle))){
        if (celloracle$GeneID[i] %in% gsgo){
            count <- count +1
            if (count == r50){
                return(i)
                break
            }
        }
    }
}

de <- read.csv("Input_X_root_markers_ROC.csv")

## Define DE
de <- read.csv("Input_X_root_markers_ROC.csv")
de <- de[c(grep("Xylem",de$cluster),grep("Phloem",de$cluster),grep("Procambium",de$cluster)),]
de <- de[!is.na(match(de$gene,exptf)),]
dat <- de %>% group_by(gene) %>% reframe(myAUC_rank = mean(myAUC_rank),avg_diff_rank = mean(avg_diff_rank))
colnames(dat) <- c("GeneID","myAUC_rank","avg_diff_rank")
de_scRNAseq <- dat

de <- read.csv("GeneSys_P_root_markers_ROC.csv")
de <- de[c(grep("Xylem",de$cluster),grep("Phloem",de$cluster),grep("Procambium",de$cluster)),]
de <- de[!is.na(match(de$gene,exptf)),]
dat <- de %>% group_by(gene) %>% reframe(myAUC_rank = mean(myAUC_rank),avg_diff_rank = mean(avg_diff_rank))
colnames(dat) <- c("GeneID","myAUC_rank","avg_diff_rank")
de_GeneSys <- dat

run_r50_de <- function(x, dat){
    de <- dat
    de$combined_score <- de[,grep(x,colnames(de))]
    de <- de %>% arrange(desc(combined_score))
    count <- 0
    for (i in seq(nrow(de))){
        if (de$GeneID[i] %in% gsgo){
            count <- count +1
            if (count == r50){
                return(i)
                break
            }
        }
    }
}

## Expressed TFs permutation
R50_permutation <- c()
for (j in 1:1000){
ran <- sample(exptf, length(exptf), replace=FALSE)
count <- 0
for (i in seq(length(ran))){
    if (ran[i] %in% gsgo){
        count <- count +1
        if (count == r50){
            R50_permutation <- c(R50_permutation,i)
            break
        }
    }
}
}

toplt <- data.frame(Methods=c("LIMA degree centrality (GeneSys)", "LIMA out centrality (GeneSys)", "LIMA in centrality (GeneSys)", "LIMA betweenness centrality (GeneSys)", 
                 "LIMA eigenvector centrality (GeneSys)","LIMA degree centrality (scRNAseq)", "LIMA out centrality (scRNAseq)", "LIMA in centrality (scRNAseq)", "LIMA betweenness centrality (scRNAseq)", 
                "LIMA eigenvector centrality (scRNAseq)", "CellOracle degree centrality (scRNAseq)", "CellOracle out centrality (scRNAseq)", 
                 "CellOracle in centrality (scRNAseq)", "CellOracle betweenness centrality (scRNAseq)", "CellOracle eigenvector centrality (scRNAseq)",
                 "CellOracle degree centrality (GeneSys)", "CellOracle out centrality (GeneSys)", 
                 "CellOracle in centrality (GeneSys)", "CellOracle betweenness centrality (GeneSys)", "CellOracle eigenvector centrality (GeneSys)",
                 "DE myAUC rank (scRNAseq)", "DE avg diff rank (scRNAseq)", "DE myAUC rank (GeneSys)", "DE avg diff rank (GeneSys)", "Permutation"),
           R50=c(run_r50_genesys('degree_centrality'),run_r50_genesys('out_centrality'),run_r50_genesys('in_centrality'),
                 run_r50_genesys('betweenness_centrality'),run_r50_genesys('eigenvector_centrality'),
                 run_r50_scRNAseq('degree_centrality'),run_r50_scRNAseq('out_centrality'),run_r50_scRNAseq('in_centrality'),
                 run_r50_scRNAseq('betweenness_centrality'),run_r50_scRNAseq('eigenvector_centrality'),
                 run_r50_celloracle('degree_centrality', celloracle_scRNAseq),run_r50_celloracle('out_centrality', celloracle_scRNAseq),run_r50_celloracle('in_centrality', celloracle_scRNAseq),
                 run_r50_celloracle('betweenness_centrality', celloracle_scRNAseq),run_r50_celloracle('eigenvector_centrality', celloracle_scRNAseq),
                 run_r50_celloracle('degree_centrality', celloracle_GeneSys),run_r50_celloracle('out_centrality', celloracle_GeneSys),run_r50_celloracle('in_centrality', celloracle_GeneSys),
                 run_r50_celloracle('betweenness_centrality', celloracle_GeneSys),run_r50_celloracle('eigenvector_centrality', celloracle_GeneSys),
                 run_r50_de('myAUC_rank', de_scRNAseq),run_r50_de('avg_diff_rank', de_scRNAseq),run_r50_de('myAUC_rank', de_GeneSys),run_r50_de('avg_diff_rank', de_GeneSys), mean(R50_permutation)))

toplt

options(repr.plot.width=10, repr.plot.height=8)
ggplot(toplt, aes(x=reorder(Methods, R50, decreasing = TRUE),  y=R50)) + geom_point(size=4)+ 
  labs(title="Stele-specific TF Prioritization Performance",x="", y = "R50")+
  theme_classic()+
  theme(axis.text.x = element_text(color = "black", 
                           size = 16, angle = 0),
          axis.text.y = element_text( color = "black", 
                           size = 16, angle = 0),
        plot.title = element_text(color="black", size=18, face="bold"),
        axis.title.x = element_text(color="black", size=16, face="bold"))+ coord_flip()

toplt <- data.frame(Methods=factor(c(rep("LIMA (GeneSys)",5), rep("LIMA (scRNAseq)",5), rep("CellOracle (scRNAseq)",5), rep("CellOracle (GeneSys)",5), rep("DE (scRNAseq)",2), rep("DE (GeneSys)",2), "Permutation"), 
                                   levels=c("Permutation", "CellOracle (scRNAseq)", "DE (scRNAseq)", "LIMA (scRNAseq)", "CellOracle (GeneSys)", "DE (GeneSys)", "LIMA (GeneSys)")), R50=toplt$R50)

options(repr.plot.width=6, repr.plot.height=5)
ggplot(toplt, aes(x=Methods, y=R50))+
  geom_boxplot(width=0.4)+ 
  labs(title="Stele-specific TF Prioritization Performance",x="", y = "R50")+
  theme_classic()+
  theme(axis.text.x = element_text(color = "black", 
                           size = 16, angle = 0),
          axis.text.y = element_text( color = "black", 
                           size = 16, angle = 0),
        plot.title = element_text(color="black", size=18, face="bold"),
        axis.title.x = element_text(color="black", size=16, face="bold"))+ coord_flip()

gsgo <- read.csv("./Gold_Standard_Root_TF_StringDB.csv", header=TRUE)
gsgo <- gsgo[which(gsgo[,6]=="Yes"),]
gsgo$preferredName
gsgo <- gsgo$GeneID
gsgo <- intersect(gsgo, gsgo_miniex_stringdb)

## R35
length(gsgo)
r50 <- 13

## Define genesys
run_r50_genesys <- function(x){
    genesys$ct_score <- min_max_normalize(rowSums(apply(genesys[,grep(paste0('atri_',x,'|tri_',x,'|lrc_',x),colnames(genesys))],2,as.numeric)))
    genesys$combined_score <- genesys$ct_score
    
    genesys <- genesys %>% arrange(desc(combined_score))
    count <- 0
    for (i in seq(nrow(genesys))){
        if (genesys$GeneID[i] %in% gsgo){
            count <- count +1
            if (count == r50){
                return(i)
                break
            }
        }
    }
}

## Define raw
run_r50_scRNAseq <- function(x){
    raw$ct_score <- min_max_normalize(rowSums(apply(raw[,grep(paste0('atri_',x,'|tri_',x,'|lrc_',x),colnames(raw))],2,as.numeric)))
    raw$combined_score <- raw$ct_score
    
    raw <- raw %>% arrange(desc(combined_score))
    count <- 0
    for (i in seq(nrow(raw))){
        if (raw$GeneID[i] %in% gsgo){
            count <- count +1
            if (count == r50){
                return(i)
                break
            }
        }
    }
}

## Define celloracle
atri <- prepros("../celloracle/atrichoblast_X_root_integrated_celloracle_gene_score_iGRN.csv")
tri <- prepros("../celloracle/trichoblast_X_root_integrated_celloracle_gene_score_iGRN.csv")
lrc <- prepros("../celloracle/lrc_X_root_integrated_celloracle_gene_score_iGRN.csv")
## Define celloracle
dat <- rbind(atri, tri, lrc)
## Remove those not present in the dataset
dat <- dat[!is.na(match(dat$X,exptf)),]
dat <- dat %>% group_by(cluster) %>% reframe(GeneID=X,degree_centrality = min_max_normalize(degree_centrality_all), 
                                  in_centrality = min_max_normalize(degree_centrality_in), 
                                  out_centrality = min_max_normalize(degree_centrality_out),
                                 betweenness_centrality = min_max_normalize(betweenness_centrality),
                                 eigenvector_centrality = min_max_normalize(eigenvector_centrality))
dat <- dat %>% group_by(GeneID) %>% reframe(degree_centrality = sum(degree_centrality), 
                                  in_centrality = sum(in_centrality), 
                                  out_centrality = sum(out_centrality),
                                 betweenness_centrality = sum(betweenness_centrality),
                                 eigenvector_centrality = sum(eigenvector_centrality))
dat <- cbind(dat[,1],apply(dat[,2:6],2,min_max_normalize))

celloracle_scRNAseq <- dat

## Define celloracle
atri <- prepros("../celloracle/atrichoblast_P_root_integrated_celloracle_gene_score_iGRN.csv")
tri <- prepros("../celloracle/trichoblast_P_root_integrated_celloracle_gene_score_iGRN.csv")
lrc <- prepros("../celloracle/lrc_P_root_integrated_celloracle_gene_score_iGRN.csv")
## Define celloracle
dat <- rbind(atri, tri, lrc)
## Remove those not present in the dataset
dat <- dat[!is.na(match(dat$X,exptf)),]
dat <- dat %>% group_by(cluster) %>% reframe(GeneID=X,degree_centrality = min_max_normalize(degree_centrality_all), 
                                  in_centrality = min_max_normalize(degree_centrality_in), 
                                  out_centrality = min_max_normalize(degree_centrality_out),
                                 betweenness_centrality = min_max_normalize(betweenness_centrality),
                                 eigenvector_centrality = min_max_normalize(eigenvector_centrality))
dat <- dat %>% group_by(GeneID) %>% reframe(degree_centrality = sum(degree_centrality), 
                                  in_centrality = sum(in_centrality), 
                                  out_centrality = sum(out_centrality),
                                 betweenness_centrality = sum(betweenness_centrality),
                                 eigenvector_centrality = sum(eigenvector_centrality))
dat <- cbind(dat[,1],apply(dat[,2:6],2,min_max_normalize))

celloracle_GeneSys <- dat

run_r50_celloracle <- function(x, dat){
    celloracle <- dat
    celloracle$combined_score <- celloracle[,grep(x,colnames(celloracle))]
    celloracle <- celloracle %>% arrange(desc(combined_score))
    count <- 0
    for (i in seq(nrow(celloracle))){
        if (celloracle$GeneID[i] %in% gsgo){
            count <- count +1
            if (count == r50){
                return(i)
                break
            }
        }
    }
}

## Define DE
de <- read.csv("Input_X_root_markers_ROC.csv")
de <- de[c(grep("Atrichoblast",de$cluster),grep("Trichoblast",de$cluster),grep("Lateral Root Cap",de$cluster)),]
de <- de[!is.na(match(de$gene,exptf)),]
dat <- de %>% group_by(gene) %>% reframe(myAUC_rank = mean(myAUC_rank),avg_diff_rank = mean(avg_diff_rank))
colnames(dat) <- c("GeneID","myAUC_rank","avg_diff_rank")
de_scRNAseq <- dat

de <- read.csv("GeneSys_P_root_markers_ROC.csv")
de <- de[c(grep("Atrichoblast",de$cluster),grep("Trichoblast",de$cluster),grep("Lateral Root Cap",de$cluster)),]
de <- de[!is.na(match(de$gene,exptf)),]
dat <- de %>% group_by(gene) %>% reframe(myAUC_rank = mean(myAUC_rank),avg_diff_rank = mean(avg_diff_rank))
colnames(dat) <- c("GeneID","myAUC_rank","avg_diff_rank")
de_GeneSys <- dat

run_r50_de <- function(x, dat){
    de <- dat
    de$combined_score <- de[,grep(x,colnames(de))]
    de <- de %>% arrange(desc(combined_score))
    count <- 0
    for (i in seq(nrow(de))){
        if (de$GeneID[i] %in% gsgo){
            count <- count +1
            if (count == r50){
                return(i)
                break
            }
        }
    }
}

## Expressed TFs permutation
R50_permutation <- c()
for (j in 1:1000){
ran <- sample(exptf, length(exptf), replace=FALSE)
count <- 0
for (i in seq(length(ran))){
    if (ran[i] %in% gsgo){
        count <- count +1
        if (count == r50){
            R50_permutation <- c(R50_permutation,i)
            break
        }
    }
}
}

length(intersect(de_GeneSys$GeneID, gsgo))

toplt <- data.frame(Methods=c("LIMA degree centrality (GeneSys)", "LIMA out centrality (GeneSys)", "LIMA in centrality (GeneSys)", "LIMA betweenness centrality (GeneSys)", 
                 "LIMA eigenvector centrality (GeneSys)","LIMA degree centrality (scRNAseq)", "LIMA out centrality (scRNAseq)", "LIMA in centrality (scRNAseq)", "LIMA betweenness centrality (scRNAseq)", 
                "LIMA eigenvector centrality (scRNAseq)", "CellOracle degree centrality (scRNAseq)", "CellOracle out centrality (scRNAseq)", 
                 "CellOracle in centrality (scRNAseq)", "CellOracle betweenness centrality (scRNAseq)", "CellOracle eigenvector centrality (scRNAseq)",
                 "CellOracle degree centrality (GeneSys)", "CellOracle out centrality (GeneSys)", 
                 "CellOracle in centrality (GeneSys)", "CellOracle betweenness centrality (GeneSys)", "CellOracle eigenvector centrality (GeneSys)",
                 "DE myAUC rank (scRNAseq)", "DE avg diff rank (scRNAseq)", "DE myAUC rank (GeneSys)", "DE avg diff rank (GeneSys)", "Permutation"),
           R50=c(run_r50_genesys('degree_centrality'),run_r50_genesys('out_centrality'),run_r50_genesys('in_centrality'),
                 run_r50_genesys('betweenness_centrality'),run_r50_genesys('eigenvector_centrality'),
                 run_r50_scRNAseq('degree_centrality'),run_r50_scRNAseq('out_centrality'),run_r50_scRNAseq('in_centrality'),
                 run_r50_scRNAseq('betweenness_centrality'),run_r50_scRNAseq('eigenvector_centrality'),
                 run_r50_celloracle('degree_centrality', celloracle_scRNAseq),run_r50_celloracle('out_centrality', celloracle_scRNAseq),run_r50_celloracle('in_centrality', celloracle_scRNAseq),
                 run_r50_celloracle('betweenness_centrality', celloracle_scRNAseq),run_r50_celloracle('eigenvector_centrality', celloracle_scRNAseq),
                 run_r50_celloracle('degree_centrality', celloracle_GeneSys),run_r50_celloracle('out_centrality', celloracle_GeneSys),run_r50_celloracle('in_centrality', celloracle_GeneSys),
                 run_r50_celloracle('betweenness_centrality', celloracle_GeneSys),run_r50_celloracle('eigenvector_centrality', celloracle_GeneSys),
                 run_r50_de('myAUC_rank', de_scRNAseq),run_r50_de('avg_diff_rank', de_scRNAseq),run_r50_de('myAUC_rank', de_GeneSys),run_r50_de('avg_diff_rank', de_GeneSys), mean(R50_permutation)))

toplt

options(repr.plot.width=10, repr.plot.height=8)
ggplot(toplt, aes(x=reorder(Methods, R50, decreasing = TRUE),  y=R50)) + geom_point(size=4)+ 
  labs(title="Epidermis-specific TF Prioritization Performance",x="", y = "R35")+
  theme_classic()+
  theme(axis.text.x = element_text(color = "black", 
                           size = 16, angle = 0),
          axis.text.y = element_text( color = "black", 
                           size = 16, angle = 0),
        plot.title = element_text(color="black", size=18, face="bold"),
        axis.title.x = element_text(color="black", size=16, face="bold"))+ coord_flip()

toplt <- data.frame(Methods=factor(c(rep("LIMA (GeneSys)",5), rep("LIMA (scRNAseq)",5), rep("CellOracle (scRNAseq)",5), rep("CellOracle (GeneSys)",5), rep("DE (scRNAseq)",2), rep("DE (GeneSys)",2), "Permutation"), 
                                   levels=c("Permutation", "CellOracle (scRNAseq)", "LIMA (scRNAseq)", "DE (scRNAseq)", "CellOracle (GeneSys)", "LIMA (GeneSys)", "DE (GeneSys)")), R50=toplt$R50)

options(repr.plot.width=6, repr.plot.height=5)
ggplot(toplt, aes(x=Methods, y=R50))+
  geom_boxplot(width=0.4)+ 
  labs(title="Epidermis-specific TF Prioritization Performance",x="", y = "R35")+
  theme_classic()+
  theme(axis.text.x = element_text(color = "black", 
                           size = 16, angle = 0),
          axis.text.y = element_text( color = "black", 
                           size = 16, angle = 0),
        plot.title = element_text(color="black", size=18, face="bold"),
        axis.title.x = element_text(color="black", size=16, face="bold"))+ coord_flip()

gsgo <- read.csv("./Gold_Standard_Root_TF_StringDB.csv", header=TRUE)
gsgo <- gsgo[which(gsgo[,7]=="Yes"),]
gsgo$preferredName
gsgo <- gsgo$GeneID
gsgo <- intersect(gsgo, gsgo_miniex_stringdb)

length(gsgo)
r50 <- 10

## Define genesys
run_r50_genesys <- function(x){
    genesys$ct_score <- min_max_normalize(rowSums(apply(genesys[,grep(paste0('xyl_',x),colnames(genesys))],2,as.numeric)))
    genesys$combined_score <- genesys$ct_score
    
    genesys <- genesys %>% arrange(desc(combined_score))
    count <- 0
    for (i in seq(nrow(genesys))){
        if (genesys$GeneID[i] %in% gsgo){
            count <- count +1
            if (count == r50){
                return(i)
                break
            }
        }
    }
}

## Define genesys
run_r50_scRNAseq <- function(x){
    raw$ct_score <- min_max_normalize(rowSums(apply(raw[,grep(paste0('xyl_',x),colnames(raw))],2,as.numeric)))
    raw$combined_score <- raw$ct_score
    
    raw <- raw %>% arrange(desc(combined_score))
    count <- 0
    for (i in seq(nrow(raw))){
        if (raw$GeneID[i] %in% gsgo){
            count <- count +1
            if (count == r50){
                return(i)
                break
            }
        }
    }
}

## Define celloracle
xyl <- prepros("../celloracle/xylem_X_root_integrated_celloracle_gene_score_iGRN.csv")
## Define celloracle
dat <- rbind(xyl)
## Remove those not present in the dataset
dat <- dat[!is.na(match(dat$X,exptf)),]
dat <- dat %>% group_by(cluster) %>% reframe(GeneID=X,degree_centrality = min_max_normalize(degree_centrality_all), 
                                  in_centrality = min_max_normalize(degree_centrality_in), 
                                  out_centrality = min_max_normalize(degree_centrality_out),
                                 betweenness_centrality = min_max_normalize(betweenness_centrality),
                                 eigenvector_centrality = min_max_normalize(eigenvector_centrality))
dat <- dat %>% group_by(GeneID) %>% reframe(degree_centrality = sum(degree_centrality), 
                                  in_centrality = sum(in_centrality), 
                                  out_centrality = sum(out_centrality),
                                 betweenness_centrality = sum(betweenness_centrality),
                                 eigenvector_centrality = sum(eigenvector_centrality))
dat <- cbind(dat[,1],apply(dat[,2:6],2,min_max_normalize))

celloracle_scRNAseq <- dat

## Define celloracle
xyl <- prepros("../celloracle/xylem_P_root_integrated_celloracle_gene_score_iGRN.csv")
## Define celloracle
dat <- rbind(xyl)
## Remove those not present in the dataset
dat <- dat[!is.na(match(dat$X,exptf)),]
dat <- dat %>% group_by(cluster) %>% reframe(GeneID=X,degree_centrality = min_max_normalize(degree_centrality_all), 
                                  in_centrality = min_max_normalize(degree_centrality_in), 
                                  out_centrality = min_max_normalize(degree_centrality_out),
                                 betweenness_centrality = min_max_normalize(betweenness_centrality),
                                 eigenvector_centrality = min_max_normalize(eigenvector_centrality))
dat <- dat %>% group_by(GeneID) %>% reframe(degree_centrality = sum(degree_centrality), 
                                  in_centrality = sum(in_centrality), 
                                  out_centrality = sum(out_centrality),
                                 betweenness_centrality = sum(betweenness_centrality),
                                 eigenvector_centrality = sum(eigenvector_centrality))
dat <- cbind(dat[,1],apply(dat[,2:6],2,min_max_normalize))

celloracle_GeneSys <- dat

run_r50_celloracle <- function(x, dat){
    celloracle <- dat
    celloracle$combined_score <- celloracle[,grep(x,colnames(celloracle))]
    celloracle <- celloracle %>% arrange(desc(combined_score))
    count <- 0
    for (i in seq(nrow(celloracle))){
        if (celloracle$GeneID[i] %in% gsgo){
            count <- count +1
            if (count == r50){
                return(i)
                break
            }
        }
    }
}

## Define DE
de <- read.csv("Input_X_root_markers_ROC.csv")
de <- de[c(grep("Xylem",de$cluster)),]
de <- de[!is.na(match(de$gene,exptf)),]
dat <- de %>% group_by(gene) %>% reframe(myAUC_rank = mean(myAUC_rank),avg_diff_rank = mean(avg_diff_rank))
colnames(dat) <- c("GeneID","myAUC_rank","avg_diff_rank")
de_scRNAseq <- dat

de <- read.csv("GeneSys_P_root_markers_ROC.csv")
de <- de[c(grep("Xylem",de$cluster)),]
de <- de[!is.na(match(de$gene,exptf)),]
dat <- de %>% group_by(gene) %>% reframe(myAUC_rank = mean(myAUC_rank),avg_diff_rank = mean(avg_diff_rank))
colnames(dat) <- c("GeneID","myAUC_rank","avg_diff_rank")
de_GeneSys <- dat

run_r50_de <- function(x, dat){
    de <- dat
    de$combined_score <- de[,grep(x,colnames(de))]
    de <- de %>% arrange(desc(combined_score))
    count <- 0
    for (i in seq(nrow(de))){
        if (de$GeneID[i] %in% gsgo){
            count <- count +1
            if (count == r50){
                return(i)
                break
            }
        }
    }
}

## Expressed TFs permutation
R50_permutation <- c()
for (j in 1:1000){
ran <- sample(exptf, length(exptf), replace=FALSE)
count <- 0
for (i in seq(length(ran))){
    if (ran[i] %in% gsgo){
        count <- count +1
        if (count == r50){
            R50_permutation <- c(R50_permutation,i)
            break
        }
    }
}
}

toplt <- data.frame(Methods=c("LIMA degree centrality (GeneSys)", "LIMA out centrality (GeneSys)", "LIMA in centrality (GeneSys)", "LIMA betweenness centrality (GeneSys)", 
                 "LIMA eigenvector centrality (GeneSys)","LIMA degree centrality (scRNAseq)", "LIMA out centrality (scRNAseq)", "LIMA in centrality (scRNAseq)", "LIMA betweenness centrality (scRNAseq)", 
                "LIMA eigenvector centrality (scRNAseq)", "CellOracle degree centrality (scRNAseq)", "CellOracle out centrality (scRNAseq)", 
                 "CellOracle in centrality (scRNAseq)", "CellOracle betweenness centrality (scRNAseq)", "CellOracle eigenvector centrality (scRNAseq)",
                 "CellOracle degree centrality (GeneSys)", "CellOracle out centrality (GeneSys)", 
                 "CellOracle in centrality (GeneSys)", "CellOracle betweenness centrality (GeneSys)", "CellOracle eigenvector centrality (GeneSys)",
                 "DE myAUC rank (scRNAseq)", "DE avg diff rank (scRNAseq)", "DE myAUC rank (GeneSys)", "DE avg diff rank (GeneSys)", "Permutation"),
           R50=c(run_r50_genesys('degree_centrality'),run_r50_genesys('out_centrality'),run_r50_genesys('in_centrality'),
                 run_r50_genesys('betweenness_centrality'),run_r50_genesys('eigenvector_centrality'),
                 run_r50_scRNAseq('degree_centrality'),run_r50_scRNAseq('out_centrality'),run_r50_scRNAseq('in_centrality'),
                 run_r50_scRNAseq('betweenness_centrality'),run_r50_scRNAseq('eigenvector_centrality'),
                 run_r50_celloracle('degree_centrality', celloracle_scRNAseq),run_r50_celloracle('out_centrality', celloracle_scRNAseq),run_r50_celloracle('in_centrality', celloracle_scRNAseq),
                 run_r50_celloracle('betweenness_centrality', celloracle_scRNAseq),run_r50_celloracle('eigenvector_centrality', celloracle_scRNAseq),
                 run_r50_celloracle('degree_centrality', celloracle_GeneSys),run_r50_celloracle('out_centrality', celloracle_GeneSys),run_r50_celloracle('in_centrality', celloracle_GeneSys),
                 run_r50_celloracle('betweenness_centrality', celloracle_GeneSys),run_r50_celloracle('eigenvector_centrality', celloracle_GeneSys),
                 run_r50_de('myAUC_rank', de_scRNAseq),run_r50_de('avg_diff_rank', de_scRNAseq),run_r50_de('myAUC_rank', de_GeneSys),run_r50_de('avg_diff_rank', de_GeneSys), mean(R50_permutation)))

toplt

options(repr.plot.width=10, repr.plot.height=8)
ggplot(toplt, aes(x=reorder(Methods, R50, decreasing = TRUE),  y=R50)) + geom_point(size=4)+ 
  labs(title="Xylem-specific TF Prioritization Performance",x="", y = "R50")+
  theme_classic()+
  theme(axis.text.x = element_text(color = "black", 
                           size = 16, angle = 0),
          axis.text.y = element_text( color = "black", 
                           size = 16, angle = 0),
        plot.title = element_text(color="black", size=18, face="bold"),
        axis.title.x = element_text(color="black", size=16, face="bold"))+ coord_flip()

toplt <- data.frame(Methods=factor(c(rep("LIMA (GeneSys)",5), rep("LIMA (scRNAseq)",5), rep("CellOracle (scRNAseq)",5), rep("CellOracle (GeneSys)",5), rep("DE (scRNAseq)",2), rep("DE (GeneSys)",2), "Permutation"), 
                                   levels=c("Permutation", "CellOracle (scRNAseq)", "DE (scRNAseq)", "LIMA (scRNAseq)", "CellOracle (GeneSys)", "DE (GeneSys)", "LIMA (GeneSys)")), R50=toplt$R50)

options(repr.plot.width=6, repr.plot.height=5)
ggplot(toplt, aes(x=Methods, y=R50))+
  geom_boxplot(width=0.4)+ 
  labs(title="Xylem-specific TF Prioritization Performance",x="", y = "R50")+
  theme_classic()+
  theme(axis.text.x = element_text(color = "black", 
                           size = 16, angle = 0),
          axis.text.y = element_text( color = "black", 
                           size = 16, angle = 0),
        plot.title = element_text(color="black", size=18, face="bold"),
        axis.title.x = element_text(color="black", size=16, face="bold"))+ coord_flip()

gsgo <- read.csv("./Gold_Standard_Root_TF_StringDB.csv", header=TRUE)
gsgo <- gsgo[which(gsgo[,8]=="Yes"),]
gsgo$preferredName
gsgo <- gsgo$GeneID
gsgo <- intersect(gsgo, gsgo_miniex_stringdb)

length(gsgo)
r50 <- 5

## Define genesys
run_r50_genesys <- function(x){
    genesys$ct_score <- min_max_normalize(rowSums(apply(genesys[,grep(paste0('^tri_',x),colnames(genesys))],2,as.numeric)))
    genesys$combined_score <- genesys$ct_score
    
    genesys <- genesys %>% arrange(desc(combined_score))
    count <- 0
    for (i in seq(nrow(genesys))){
        if (genesys$GeneID[i] %in% gsgo){
            count <- count +1
            if (count == r50){
                return(i)
                break
            }
        }
    }
}

## Define raw
run_r50_scRNAseq <- function(x){
    raw$ct_score <- min_max_normalize(rowSums(apply(raw[,grep(paste0('^tri_',x),colnames(raw))],2,as.numeric)))
    raw$combined_score <- raw$ct_score
    
    raw <- raw %>% arrange(desc(combined_score))
    count <- 0
    for (i in seq(nrow(raw))){
        if (raw$GeneID[i] %in% gsgo){
            count <- count +1
            if (count == r50){
                return(i)
                break
            }
        }
    }
}

## Define celloracle
tri <- prepros("../celloracle/trichoblast_X_root_integrated_celloracle_gene_score_iGRN.csv")
## Define celloracle
dat <- rbind(tri)
## Remove those not present in the dataset
dat <- dat[!is.na(match(dat$X,exptf)),]
dat <- dat %>% group_by(cluster) %>% reframe(GeneID=X,degree_centrality = min_max_normalize(degree_centrality_all), 
                                  in_centrality = min_max_normalize(degree_centrality_in), 
                                  out_centrality = min_max_normalize(degree_centrality_out),
                                 betweenness_centrality = min_max_normalize(betweenness_centrality),
                                 eigenvector_centrality = min_max_normalize(eigenvector_centrality))
dat <- dat %>% group_by(GeneID) %>% reframe(degree_centrality = sum(degree_centrality), 
                                  in_centrality = sum(in_centrality), 
                                  out_centrality = sum(out_centrality),
                                 betweenness_centrality = sum(betweenness_centrality),
                                 eigenvector_centrality = sum(eigenvector_centrality))
dat <- cbind(dat[,1],apply(dat[,2:6],2,min_max_normalize))

celloracle_scRNAseq <- dat

## Define celloracle
tri <- prepros("../celloracle/trichoblast_P_root_integrated_celloracle_gene_score_iGRN.csv")
## Define celloracle
dat <- rbind(tri)
## Remove those not present in the dataset
dat <- dat[!is.na(match(dat$X,exptf)),]
dat <- dat %>% group_by(cluster) %>% reframe(GeneID=X,degree_centrality = min_max_normalize(degree_centrality_all), 
                                  in_centrality = min_max_normalize(degree_centrality_in), 
                                  out_centrality = min_max_normalize(degree_centrality_out),
                                 betweenness_centrality = min_max_normalize(betweenness_centrality),
                                 eigenvector_centrality = min_max_normalize(eigenvector_centrality))
dat <- dat %>% group_by(GeneID) %>% reframe(degree_centrality = sum(degree_centrality), 
                                  in_centrality = sum(in_centrality), 
                                  out_centrality = sum(out_centrality),
                                 betweenness_centrality = sum(betweenness_centrality),
                                 eigenvector_centrality = sum(eigenvector_centrality))
dat <- cbind(dat[,1],apply(dat[,2:6],2,min_max_normalize))

celloracle_GeneSys <- dat

run_r50_celloracle <- function(x, dat){
    celloracle <- dat
    celloracle$combined_score <- celloracle[,grep(x,colnames(celloracle))]
    celloracle <- celloracle %>% arrange(desc(combined_score))
    count <- 0
    for (i in seq(nrow(celloracle))){
        if (celloracle$GeneID[i] %in% gsgo){
            count <- count +1
            if (count == r50){
                return(i)
                break
            }
        }
    }
}

## Define DE
de <- read.csv("Input_X_root_markers_ROC.csv")
de <- de[c(grep("Trichoblast",de$cluster)),]
de <- de[!is.na(match(de$gene,exptf)),]
dat <- de %>% group_by(gene) %>% reframe(myAUC_rank = mean(myAUC_rank),avg_diff_rank = mean(avg_diff_rank))
colnames(dat) <- c("GeneID","myAUC_rank","avg_diff_rank")
de_scRNAseq <- dat

de <- read.csv("GeneSys_P_root_markers_ROC.csv")
de <- de[c(grep("Trichoblast",de$cluster)),]
de <- de[!is.na(match(de$gene,exptf)),]
dat <- de %>% group_by(gene) %>% reframe(myAUC_rank = mean(myAUC_rank),avg_diff_rank = mean(avg_diff_rank))
colnames(dat) <- c("GeneID","myAUC_rank","avg_diff_rank")
de_GeneSys <- dat

run_r50_de <- function(x, dat){
    de <- dat
    de$combined_score <- de[,grep(x,colnames(de))]
    de <- de %>% arrange(desc(combined_score))
    count <- 0
    for (i in seq(nrow(de))){
        if (de$GeneID[i] %in% gsgo){
            count <- count +1
            if (count == r50){
                return(i)
                break
            }
        }
    }
}

length(intersect(de_GeneSys$GeneID, gsgo))

## Expressed TFs permutation
R50_permutation <- c()
for (j in 1:1000){
ran <- sample(exptf, length(exptf), replace=FALSE)
count <- 0
for (i in seq(length(ran))){
    if (ran[i] %in% gsgo){
        count <- count +1
        if (count == r50){
            R50_permutation <- c(R50_permutation,i)
            break
        }
    }
}
}

toplt <- data.frame(Methods=c("LIMA degree centrality (GeneSys)", "LIMA out centrality (GeneSys)", "LIMA in centrality (GeneSys)", "LIMA betweenness centrality (GeneSys)", 
                 "LIMA eigenvector centrality (GeneSys)","LIMA degree centrality (scRNAseq)", "LIMA out centrality (scRNAseq)", "LIMA in centrality (scRNAseq)", "LIMA betweenness centrality (scRNAseq)", 
                "LIMA eigenvector centrality (scRNAseq)", "CellOracle degree centrality (scRNAseq)", "CellOracle out centrality (scRNAseq)", 
                 "CellOracle in centrality (scRNAseq)", "CellOracle betweenness centrality (scRNAseq)", "CellOracle eigenvector centrality (scRNAseq)",
                 "CellOracle degree centrality (GeneSys)", "CellOracle out centrality (GeneSys)", 
                 "CellOracle in centrality (GeneSys)", "CellOracle betweenness centrality (GeneSys)", "CellOracle eigenvector centrality (GeneSys)",
                 "DE myAUC rank (scRNAseq)", "DE avg diff rank (scRNAseq)", "DE myAUC rank (GeneSys)", "DE avg diff rank (GeneSys)", "Permutation"),
           R50=c(run_r50_genesys('degree_centrality'),run_r50_genesys('out_centrality'),run_r50_genesys('in_centrality'),
                 run_r50_genesys('betweenness_centrality'),run_r50_genesys('eigenvector_centrality'),
                 run_r50_scRNAseq('degree_centrality'),run_r50_scRNAseq('out_centrality'),run_r50_scRNAseq('in_centrality'),
                 run_r50_scRNAseq('betweenness_centrality'),run_r50_scRNAseq('eigenvector_centrality'),
                 run_r50_celloracle('degree_centrality', celloracle_scRNAseq),run_r50_celloracle('out_centrality', celloracle_scRNAseq),run_r50_celloracle('in_centrality', celloracle_scRNAseq),
                 run_r50_celloracle('betweenness_centrality', celloracle_scRNAseq),run_r50_celloracle('eigenvector_centrality', celloracle_scRNAseq),
                 run_r50_celloracle('degree_centrality', celloracle_GeneSys),run_r50_celloracle('out_centrality', celloracle_GeneSys),run_r50_celloracle('in_centrality', celloracle_GeneSys),
                 run_r50_celloracle('betweenness_centrality', celloracle_GeneSys),run_r50_celloracle('eigenvector_centrality', celloracle_GeneSys),
                 run_r50_de('myAUC_rank', de_scRNAseq),run_r50_de('avg_diff_rank', de_scRNAseq),run_r50_de('myAUC_rank', de_GeneSys),run_r50_de('avg_diff_rank', de_GeneSys), mean(R50_permutation)))

toplt

options(repr.plot.width=10, repr.plot.height=8)
ggplot(toplt, aes(x=reorder(Methods, R50, decreasing = TRUE),  y=R50)) + geom_point(size=4)+ 
  labs(title="Trichoblast-specific TF Prioritization Performance",x="", y = "R15")+
  theme_classic()+
  theme(axis.text.x = element_text(color = "black", 
                           size = 16, angle = 0),
          axis.text.y = element_text( color = "black", 
                           size = 16, angle = 0),
        plot.title = element_text(color="black", size=18, face="bold"),
        axis.title.x = element_text(color="black", size=16, face="bold"))+ coord_flip()

toplt <- data.frame(Methods=factor(c(rep("LIMA (GeneSys)",5), rep("LIMA (scRNAseq)",5), rep("CellOracle (scRNAseq)",5), rep("CellOracle (GeneSys)",5), rep("DE (scRNAseq)",2), rep("DE (GeneSys)",2), "Permutation"), 
                                   levels=c("Permutation", "CellOracle (scRNAseq)", "DE (scRNAseq)", "LIMA (scRNAseq)", "CellOracle (GeneSys)", "DE (GeneSys)", "LIMA (GeneSys)")), R50=toplt$R50)

options(repr.plot.width=6, repr.plot.height=5)
ggplot(toplt, aes(x=Methods, y=R50))+
  geom_boxplot(width=0.4)+ 
  labs(title="Trichoblast-specific TF Prioritization Performance",x="", y = "R15")+
  theme_classic()+
  theme(axis.text.x = element_text(color = "black", 
                           size = 16, angle = 0),
          axis.text.y = element_text( color = "black", 
                           size = 16, angle = 0),
        plot.title = element_text(color="black", size=18, face="bold"),
        axis.title.x = element_text(color="black", size=16, face="bold"))+ coord_flip()

dat <- genesys

plot_heatmap <- function(gene, centrality){
    sub <- as.numeric(t(dat[which(dat$TF==gene),grep(centrality,colnames(dat))]))
    sub <- as.matrix(cbind(sub[1:10],sub[11:20],sub[21:30],sub[31:40],sub[41:50]))
    ## Cell types
    rownames(sub) <- c('Atrichoblast','Trichoblast','LRC','Cortex','Endodermis','Pericycle','Procambium','Xylem','Phloem','Columella')
    ## Transition
    colnames(sub) <- c('t0-t1', 't1-t3', 't3-t5', 't5-t7', 't7-t9')
    # Reorder rows
    sub <- sub[c(8,9,7,6,5,4,2,1,3,10),]
    col_fun = colorRamp2(c(0, 0.001, 1), c('navy',"khaki1", "red"))

    Heatmap(sub, rect_gp = gpar(col = "white", lwd = 1) , show_heatmap_legend = TRUE, name = paste0(gene,"\n","weighted","\n","network","\n","centrality"),
            col = col_fun, column_title = paste0(str_split_i(centrality,"_",1),"\n",str_split_i(centrality,"_",2)), column_title_gp = gpar(fontsize = 12, fontface = "bold"), 
            cluster_rows = FALSE, cluster_columns = FALSE,
            show_column_names = FALSE, bottom_annotation = HeatmapAnnotation(
            text = anno_text(colnames(sub), rot = 60, location = unit(1, "npc"), just = "right"))
       )
}

plot_heatmap2 <- function(gene, centrality){
    sub <- as.numeric(t(dat[which(dat$TF==gene),grep(centrality,colnames(dat))]))
    sub <- as.matrix(cbind(sub[1:10],sub[11:20],sub[21:30],sub[31:40],sub[41:50]))
    ## Cell types
    rownames(sub) <- c('Atrichoblast','Trichoblast','LRC','Cortex','Endodermis','Pericycle','Procambium','Xylem','Phloem','Columella')
    ## Transition
    colnames(sub) <- c('t0-t1', 't1-t3', 't3-t5', 't5-t7', 't7-t9')
    # Reorder rows
    sub <- sub[c(8,9,7,6,5,4,2,1,3,10),]
    col_fun = colorRamp2(c(0, 0.001, 1), c('gainsboro', "grey", "black"))

    Heatmap(sub, rect_gp = gpar(col = "white", lwd = 1) , show_heatmap_legend = FALSE,
            col = col_fun, column_title = "out degree \n centrality", column_title_gp = gpar(fontsize = 12, fontface = "bold"), 
            cluster_rows = FALSE, cluster_columns = FALSE,
            show_column_names = FALSE, bottom_annotation = HeatmapAnnotation(
            text = anno_text(colnames(sub), rot = 60, location = unit(1, "npc"), just = "right"))
       )
}

plot_heatmap3 <- function(gene, centrality){
    sub <- as.numeric(t(dat[which(dat$TF==gene),grep(centrality,colnames(dat))]))
    sub <- as.matrix(cbind(sub[1:10],sub[11:20],sub[21:30],sub[31:40],sub[41:50]))
    ## Cell types
    rownames(sub) <- c('Atrichoblast','Trichoblast','LRC','Cortex','Endodermis','Pericycle','Procambium','Xylem','Phloem','Columella')
    ## Transition
    colnames(sub) <- c('t0-t1', 't1-t3', 't3-t5', 't5-t7', 't7-t9')
    # Reorder rows
    sub <- sub[c(8,9,7,6,5,4,2,1,3,10),]
    col_fun = colorRamp2(c(0, 0.001, 1), c('gainsboro', "grey", "black"))

    Heatmap(sub, rect_gp = gpar(col = "white", lwd = 1), name = paste0("unweighted","\n","network","\n","centrality"),
            col = col_fun, column_title = "in degree \n centrality", column_title_gp = gpar(fontsize = 12, fontface = "bold"), 
            cluster_rows = FALSE, cluster_columns = FALSE,
            show_column_names = FALSE, bottom_annotation = HeatmapAnnotation(
            text = anno_text(colnames(sub), rot = 60, location = unit(1, "npc"), just = "right"))
       )
}

plot_all_centrality <- function(gene){
    options(repr.plot.width=10, repr.plot.height=6)
    plot_heatmap(gene,"betweenness_centrality") + plot_heatmap2(gene,"out_centrality") + plot_heatmap3(gene,"in_centrality")
}

plot_all_centrality("TTG1")

plot_all_centrality("CPC")

plot_all_centrality("SHR")

plot_all_centrality("SCR")

plot_all_centrality("BLJ")

plot_all_centrality("JKD")

plot_all_centrality("MYB36")

plot_all_centrality("RVN")

plot_all_centrality("MGP")

plot_all_centrality("NUC")

plot_all_centrality("WER")

## HAT7
plot_all_centrality("HAT7")

## GATA10
plot_all_centrality("GATA10")

## GATA11
plot_all_centrality("GATA11")

plot_all_centrality("AN3")

plot_all_centrality("GL2")

plot_all_centrality("LBD15")

plot_all_centrality <- function(gene){
    options(repr.plot.width=8, repr.plot.height=4)
    plot_heatmap(gene,"betweenness_centrality") + plot_heatmap2(gene,"out_centrality") + plot_heatmap3(gene,"in_centrality")
}

plot_all_centrality("SHR")

plot_all_centrality("WER")

plot_bc <- function(gene){
    sub <- as.numeric(t(dat[which(dat$TF==gene),grep("betweenness_centrality",colnames(dat))]))
    sub <- as.matrix(cbind(sub[1:10],sub[11:20],sub[21:30],sub[31:40],sub[41:50]))
    ## Cell types
    rownames(sub) <- c('Atrichoblast','Trichoblast','LRC','Cortex','Endodermis','Pericycle','Procambium','Xylem','Phloem','Columella')
    ## Transition
    colnames(sub) <- c('t0-t1', 't1-t3', 't3-t5', 't5-t7', 't7-t9')
    # Reorder rows
    sub <- sub[c(8,9,7,6,5,4,2,1,3,10),]
    col_fun = colorRamp2(c(0, 0.001, 1), c('navy',"khaki1", "red"))

    Heatmap(sub, rect_gp = gpar(col = "white", lwd = 1) , show_heatmap_legend = FALSE,
            col = col_fun, column_title = gene, column_title_gp = gpar(fontsize = 12, fontface = "bold"), 
            cluster_rows = FALSE, cluster_columns = FALSE,
            show_column_names = FALSE, bottom_annotation = HeatmapAnnotation(
            text = anno_text(colnames(sub), rot = 60, location = unit(1, "npc"), just = "right"))
       )
}

plot_oc <- function(gene){
    sub <- as.numeric(t(dat[which(dat$TF==gene),grep("out_centrality",colnames(dat))]))
    sub <- as.matrix(cbind(sub[1:10],sub[11:20],sub[21:30],sub[31:40],sub[41:50]))
    ## Cell types
    rownames(sub) <- c('Atrichoblast','Trichoblast','LRC','Cortex','Endodermis','Pericycle','Procambium','Xylem','Phloem','Columella')
    ## Transition
    colnames(sub) <- c('t0-t1', 't1-t3', 't3-t5', 't5-t7', 't7-t9')
    # Reorder rows
    sub <- sub[c(8,9,7,6,5,4,2,1,3,10),]
    col_fun = colorRamp2(c(0, 0.001, 1), c('gainsboro', "grey", "black"))

    Heatmap(sub, rect_gp = gpar(col = "white", lwd = 1) , show_heatmap_legend = FALSE,
            col = col_fun, column_title = gene, column_title_gp = gpar(fontsize = 12, fontface = "bold"), 
            cluster_rows = FALSE, cluster_columns = FALSE,
            show_column_names = FALSE, bottom_annotation = HeatmapAnnotation(
            text = anno_text(colnames(sub), rot = 60, location = unit(1, "npc"), just = "right"))
       )
}

plot_ic <- function(gene){
    sub <- as.numeric(t(dat[which(dat$TF==gene),grep("in_centrality",colnames(dat))]))
    sub <- as.matrix(cbind(sub[1:10],sub[11:20],sub[21:30],sub[31:40],sub[41:50]))
    ## Cell types
    rownames(sub) <- c('Atrichoblast','Trichoblast','LRC','Cortex','Endodermis','Pericycle','Procambium','Xylem','Phloem','Columella')
    ## Transition
    colnames(sub) <- c('t0-t1', 't1-t3', 't3-t5', 't5-t7', 't7-t9')
    # Reorder rows
    sub <- sub[c(8,9,7,6,5,4,2,1,3,10),]
    col_fun = colorRamp2(c(0, 0.001, 1), c('gainsboro', "grey", "black"))

    Heatmap(sub, rect_gp = gpar(col = "white", lwd = 1) , show_heatmap_legend = FALSE,
            col = col_fun, column_title = gene, column_title_gp = gpar(fontsize = 12, fontface = "bold"), 
            cluster_rows = FALSE, cluster_columns = FALSE,
            show_column_names = FALSE, bottom_annotation = HeatmapAnnotation(
            text = anno_text(colnames(sub), rot = 60, location = unit(1, "npc"), just = "right"))
       )
}

plot_dc <- function(gene){
    sub <- as.numeric(t(dat[which(dat$TF==gene),grep("degree_centrality",colnames(dat))]))
    sub <- as.matrix(cbind(sub[1:10],sub[11:20],sub[21:30],sub[31:40],sub[41:50]))
    ## Cell types
    rownames(sub) <- c('Atrichoblast','Trichoblast','LRC','Cortex','Endodermis','Pericycle','Procambium','Xylem','Phloem','Columella')
    ## Transition
    colnames(sub) <- c('t0-t1', 't1-t3', 't3-t5', 't5-t7', 't7-t9')
    # Reorder rows
    sub <- sub[c(8,9,7,6,5,4,2,1,3,10),]
    col_fun = colorRamp2(c(0, 0.001, 1), c('gainsboro', "grey", "black"))

    Heatmap(sub, rect_gp = gpar(col = "white", lwd = 1) , show_heatmap_legend = FALSE,
            col = col_fun, column_title = gene, column_title_gp = gpar(fontsize = 12, fontface = "bold"), 
            cluster_rows = FALSE, cluster_columns = FALSE,
            show_column_names = FALSE, bottom_annotation = HeatmapAnnotation(
            text = anno_text(colnames(sub), rot = 60, location = unit(1, "npc"), just = "right"))
       )
}

plot_blank <- function(gene){
    sub <- as.numeric(t(dat[which(dat$TF==gene),grep("betweenness_centrality",colnames(dat))]))
    sub <- as.matrix(cbind(sub[1:10],sub[11:20],sub[21:30],sub[31:40],sub[41:50]))
    ## Cell types
    rownames(sub) <- c('Atrichoblast','Trichoblast','LRC','Cortex','Endodermis','Pericycle','Procambium','Xylem','Phloem','Columella')
    ## Transition
    colnames(sub) <- c('t0-t1', 't1-t3', 't3-t5', 't5-t7', 't7-t9')
    # Reorder rows
    sub <- sub[c(8,9,7,6,5,4,2,1,3,10),]
    col_fun = colorRamp2(c(0, 0.001, 1), c('white',"white", "white"))

    Heatmap(sub, rect_gp = gpar(col = "white", lwd = 1) , show_heatmap_legend = FALSE,
            col = col_fun, column_title_gp = gpar(fontsize = 12, fontface = "bold"), 
            cluster_rows = FALSE, cluster_columns = FALSE,
            show_column_names = FALSE, show_row_names = TRUE
       )
}

# Combine the top three centralies 
sub <- dat[,grep("betweenness_centrality|in_centrality|out_centrality",colnames(dat))]
sub <- as.data.frame(sapply(sub, as.numeric))
rownames(sub) <- dat$TF

head(sub)

bc_rank <- data.frame(all=rowSums(sub),atri=rowSums(sub[,grep("^atri_",colnames(sub))]),tri=rowSums(sub[,grep("^tri_",colnames(sub))])
                     ,cor=rowSums(sub[,grep("^cor_",colnames(sub))]),end=rowSums(sub[,grep("^end_",colnames(sub))])
                     ,per=rowSums(sub[,grep("^per_",colnames(sub))]),pro=rowSums(sub[,grep("^pro_",colnames(sub))])
                     ,xyl=rowSums(sub[,grep("^xyl_",colnames(sub))]),phl=rowSums(sub[,grep("^phl_",colnames(sub))])
                     ,lrc=rowSums(sub[,grep("^lrc_",colnames(sub))]),col=rowSums(sub[,grep("^col_",colnames(sub))]))

bc_rank$GeneID <- wanted_TFs$GeneID[match(rownames(bc_rank),wanted_TFs$Name)]

head(bc_rank)

atri_rank <- bc_rank[which(bc_rank$atri*2 > bc_rank$all),]%>% arrange(desc(atri))
atri_rank$GeneName <- rownames(atri_rank)

atri_rank

options(repr.plot.width=6, repr.plot.height=4)
ggplot(atri_rank[1:10,], aes(x=reorder(GeneName, atri, decreasing = FALSE),  y=atri)) + geom_point(size=4)+ 
  labs(title="Atrichoblast-specific TF Prioritization",x="", y = "Combined centrality score")+
  theme_classic()+
  theme(axis.text.x = element_text(color = "black", 
                           size = 16, angle = 0),
          axis.text.y = element_text( color = "black", 
                           size = 16, angle = 0),
        plot.title = element_text(color="black", size=18, face="bold"),
        axis.title.x = element_text(color="black", size=16, face="bold"))+ coord_flip()

write.csv(atri_rank,"Atrichoblast_TF_centrality_all_transition_zscore3.csv", quote=FALSE, row.names=TRUE)

tf_rank <- atri_rank %>% rownames(.)

p1 <- plot_bc(tf_rank[1]) + plot_bc(tf_rank[2]) + plot_bc(tf_rank[3]) + plot_bc(tf_rank[4]) + plot_bc(tf_rank[5])
p2 <- plot_bc(tf_rank[6]) + plot_bc(tf_rank[7]) + plot_bc(tf_rank[8]) + plot_bc(tf_rank[9]) + plot_bc(tf_rank[10])
p3 <- plot_bc(tf_rank[11]) + plot_bc(tf_rank[12]) + plot_bc(tf_rank[13]) + plot_bc(tf_rank[14]) + plot_bc(tf_rank[15])
p4 <- plot_bc(tf_rank[16]) + plot_bc(tf_rank[17]) + plot_bc(tf_rank[18]) + plot_bc(tf_rank[19]) + plot_bc(tf_rank[20])

p5 <- plot_oc(tf_rank[1]) + plot_oc(tf_rank[2]) + plot_oc(tf_rank[3]) + plot_oc(tf_rank[4]) + plot_oc(tf_rank[5])
p6 <- plot_oc(tf_rank[6]) + plot_oc(tf_rank[7]) + plot_oc(tf_rank[8]) + plot_oc(tf_rank[9]) + plot_oc(tf_rank[10])
p7 <- plot_oc(tf_rank[11]) + plot_oc(tf_rank[12]) + plot_oc(tf_rank[13]) + plot_oc(tf_rank[14]) + plot_oc(tf_rank[15])
p8 <- plot_oc(tf_rank[16]) + plot_oc(tf_rank[17]) + plot_oc(tf_rank[18]) + plot_oc(tf_rank[19]) + plot_oc(tf_rank[20])

p9 <- plot_ic(tf_rank[1]) + plot_ic(tf_rank[2]) + plot_ic(tf_rank[3]) + plot_ic(tf_rank[4]) + plot_ic(tf_rank[5])
p10 <- plot_ic(tf_rank[6]) + plot_ic(tf_rank[7]) + plot_ic(tf_rank[8]) + plot_ic(tf_rank[9]) + plot_ic(tf_rank[10])
p11 <- plot_ic(tf_rank[11]) + plot_ic(tf_rank[12]) + plot_ic(tf_rank[13]) + plot_ic(tf_rank[14]) + plot_ic(tf_rank[15])
p12 <- plot_ic(tf_rank[16]) + plot_ic(tf_rank[17]) + plot_ic(tf_rank[18]) + plot_ic(tf_rank[19]) + plot_ic(tf_rank[20])

p1 <- grid.grabExpr(draw(p1))
p2 <- grid.grabExpr(draw(p2))
p3 <- grid.grabExpr(draw(p3))
p4 <- grid.grabExpr(draw(p4))
p5 <- grid.grabExpr(draw(p5))
p6 <- grid.grabExpr(draw(p6))
p7 <- grid.grabExpr(draw(p7))
p8 <- grid.grabExpr(draw(p8))
p9 <- grid.grabExpr(draw(p9))
p10 <- grid.grabExpr(draw(p10))
p11 <- grid.grabExpr(draw(p11))
p12 <- grid.grabExpr(draw(p12))
col_fun = colorRamp2(c(0, 0.001, 1), c('navy',"khaki1", "red"))
legend_grob1 <- grid.grabExpr(draw(Legend(col_fun = col_fun, title = "weighted\nnetwork\ncentrality")))
col_fun = colorRamp2(c(0, 0.001, 1), c('gainsboro', "grey", "black"))
legend_grob2 <- grid.grabExpr(draw(Legend(col_fun = col_fun, title = "unweighted\nnetwork\ncentrality")))

options(repr.plot.width=24, repr.plot.height=12)
plot_grid(plot_grid(textbox_grob("Main Regulatory TFs for Atrichoblast ", gp = gpar(fontsize = 20, col='black'), background_gp = gpar(fill = "#FFFFFF", col = NA)),
          plot_grid(plot_grid(textbox_grob("Betweenness centrality", gp = gpar(fontsize = 16, col='black'), background_gp = gpar(fill = "#FFFFFF", col = NA)),plot_grid(p1,p2,p3,p4,nrow=4),nrow=2,rel_heights = c(1, 14))
                    ,plot_grid(textbox_grob("Out-degree centrality", gp = gpar(fontsize = 16, col='black'), background_gp = gpar(fill = "#FFFFFF", col = NA)),plot_grid(p5,p6,p7,p8,nrow=4),nrow=2,rel_heights=c(1,14))
                    ,plot_grid(textbox_grob("In-degree centrality", gp = gpar(fontsize = 16, col='black'), background_gp = gpar(fill = "#FFFFFF", col = NA)),plot_grid(p9,p10,p11,p12,nrow=4),nrow=2,rel_heights=c(1,14)),nrow=1),nrow=2,rel_heights = c(1, 14)),plot_grid(legend_grob1,legend_grob2,nrow=2), rel_widths = c(19, 1))

q1 <- plot_bc(tf_rank[1]) + plot_bc(tf_rank[2]) + plot_bc(tf_rank[3])
q2 <- plot_oc(tf_rank[1]) + plot_oc(tf_rank[2]) + plot_oc(tf_rank[3])
q3 <- plot_ic(tf_rank[1]) + plot_ic(tf_rank[2]) + plot_ic(tf_rank[3])

q1 <- grid.grabExpr(draw(q1))
q2 <- grid.grabExpr(draw(q2))
q3 <- grid.grabExpr(draw(q3))

options(repr.plot.width=6, repr.plot.height=8)
plot_grid(plot_grid(textbox_grob("Betweenness centrality", gp = gpar(fontsize = 16, col='black'), background_gp = gpar(fill = "#FFFFFF", col = NA)),q1,nrow=2,rel_heights = c(1, 8))
          ,plot_grid(textbox_grob("Out-degree centrality", gp = gpar(fontsize = 16, col='black'), background_gp = gpar(fill = "#FFFFFF", col = NA)),q2,nrow=2,rel_heights = c(1, 8))
          ,plot_grid(textbox_grob("In-degree centrality", gp = gpar(fontsize = 16, col='black'), background_gp = gpar(fill = "#FFFFFF", col = NA)),q3,nrow=2,rel_heights = c(1, 8)),nrow=3)

tri_rank <- bc_rank[which(bc_rank$tri*2 > bc_rank$all),]%>% arrange(desc(tri))
tri_rank$GeneName <- rownames(tri_rank)

tri_rank

options(repr.plot.width=6, repr.plot.height=4)
ggplot(tri_rank[1:10,], aes(x=reorder(GeneName, tri, decreasing = FALSE),  y=tri)) + geom_point(size=4)+ 
  labs(title="Trichoblast-specific TF Prioritization",x="", y = "Combined centrality score")+
  theme_classic()+
  theme(axis.text.x = element_text(color = "black", 
                           size = 16, angle = 0),
          axis.text.y = element_text( color = "black", 
                           size = 16, angle = 0),
        plot.title = element_text(color="black", size=18, face="bold"),
        axis.title.x = element_text(color="black", size=16, face="bold"))+ coord_flip()

write.csv(tri_rank,"Trichoblast_TF_centrality_all_transition_zscore3.csv", quote=FALSE, row.names=TRUE)

## Top20 only
tf_rank <- tri_rank %>% rownames(.)
p1 <- plot_bc(tf_rank[1]) + plot_bc(tf_rank[2]) + plot_bc(tf_rank[3]) + plot_bc(tf_rank[4]) + plot_bc(tf_rank[5])
p2 <- plot_bc(tf_rank[6]) + plot_bc(tf_rank[7]) + plot_bc(tf_rank[8]) + plot_bc(tf_rank[9]) + plot_bc(tf_rank[10])
p3 <- plot_bc(tf_rank[11]) + plot_bc(tf_rank[12]) + plot_bc(tf_rank[13]) + plot_bc(tf_rank[14]) + plot_bc(tf_rank[15])
p4 <- plot_bc(tf_rank[16]) + plot_bc(tf_rank[17]) + plot_bc(tf_rank[18]) + plot_bc(tf_rank[19]) + plot_bc(tf_rank[20])

p5 <- plot_oc(tf_rank[1]) + plot_oc(tf_rank[2]) + plot_oc(tf_rank[3]) + plot_oc(tf_rank[4]) + plot_oc(tf_rank[5])
p6 <- plot_oc(tf_rank[6]) + plot_oc(tf_rank[7]) + plot_oc(tf_rank[8]) + plot_oc(tf_rank[9]) + plot_oc(tf_rank[10])
p7 <- plot_oc(tf_rank[11]) + plot_oc(tf_rank[12]) + plot_oc(tf_rank[13]) + plot_oc(tf_rank[14]) + plot_oc(tf_rank[15])
p8 <- plot_oc(tf_rank[16]) + plot_oc(tf_rank[17]) + plot_oc(tf_rank[18]) + plot_oc(tf_rank[19]) + plot_oc(tf_rank[20])

p9 <- plot_ic(tf_rank[1]) + plot_ic(tf_rank[2]) + plot_ic(tf_rank[3]) + plot_ic(tf_rank[4]) + plot_ic(tf_rank[5])
p10 <- plot_ic(tf_rank[6]) + plot_ic(tf_rank[7]) + plot_ic(tf_rank[8]) + plot_ic(tf_rank[9]) + plot_ic(tf_rank[10])
p11 <- plot_ic(tf_rank[11]) + plot_ic(tf_rank[12]) + plot_ic(tf_rank[13]) + plot_ic(tf_rank[14]) + plot_ic(tf_rank[15])
p12 <- plot_ic(tf_rank[16]) + plot_ic(tf_rank[17]) + plot_ic(tf_rank[18]) + plot_ic(tf_rank[19]) + plot_ic(tf_rank[20])

p1 <- grid.grabExpr(draw(p1))
p2 <- grid.grabExpr(draw(p2))
p3 <- grid.grabExpr(draw(p3))
p4 <- grid.grabExpr(draw(p4))
p5 <- grid.grabExpr(draw(p5))
p6 <- grid.grabExpr(draw(p6))
p7 <- grid.grabExpr(draw(p7))
p8 <- grid.grabExpr(draw(p8))
p9 <- grid.grabExpr(draw(p9))
p10 <- grid.grabExpr(draw(p10))
p11 <- grid.grabExpr(draw(p11))
p12 <- grid.grabExpr(draw(p12))
col_fun = colorRamp2(c(0, 0.001, 1), c('navy',"khaki1", "red"))
legend_grob1 <- grid.grabExpr(draw(Legend(col_fun = col_fun, title = "weighted\nnetwork\ncentrality")))
col_fun = colorRamp2(c(0, 0.001, 1), c('gainsboro', "grey", "black"))
legend_grob2 <- grid.grabExpr(draw(Legend(col_fun = col_fun, title = "unweighted\nnetwork\ncentrality")))

options(repr.plot.width=24, repr.plot.height=12)
plot_grid(plot_grid(textbox_grob("Main Regulatory TFs for Trichoblast ", gp = gpar(fontsize = 20, col='black'), background_gp = gpar(fill = "#FFFFFF", col = NA)),
          plot_grid(plot_grid(textbox_grob("Betweenness centrality", gp = gpar(fontsize = 16, col='black'), background_gp = gpar(fill = "#FFFFFF", col = NA)),plot_grid(p1,p2,p3,p4,nrow=4),nrow=2,rel_heights = c(1, 14))
                    ,plot_grid(textbox_grob("Out-degree centrality", gp = gpar(fontsize = 16, col='black'), background_gp = gpar(fill = "#FFFFFF", col = NA)),plot_grid(p5,p6,p7,p8,nrow=4),nrow=2,rel_heights=c(1,14))
                    ,plot_grid(textbox_grob("In-degree centrality", gp = gpar(fontsize = 16, col='black'), background_gp = gpar(fill = "#FFFFFF", col = NA)),plot_grid(p9,p10,p11,p12,nrow=4),nrow=2,rel_heights=c(1,14)),nrow=1),nrow=2,rel_heights = c(1, 14)),plot_grid(legend_grob1,legend_grob2,nrow=2), rel_widths = c(19, 1))

q1 <- plot_bc(tf_rank[1]) + plot_bc(tf_rank[2]) + plot_bc(tf_rank[3])
q2 <- plot_oc(tf_rank[1]) + plot_oc(tf_rank[2]) + plot_oc(tf_rank[3])
q3 <- plot_ic(tf_rank[1]) + plot_ic(tf_rank[2]) + plot_ic(tf_rank[3])

q1 <- grid.grabExpr(draw(q1))
q2 <- grid.grabExpr(draw(q2))
q3 <- grid.grabExpr(draw(q3))

options(repr.plot.width=6, repr.plot.height=8)
plot_grid(plot_grid(textbox_grob("Betweenness centrality", gp = gpar(fontsize = 16, col='black'), background_gp = gpar(fill = "#FFFFFF", col = NA)),q1,nrow=2,rel_heights = c(1, 8))
          ,plot_grid(textbox_grob("Out-degree centrality", gp = gpar(fontsize = 16, col='black'), background_gp = gpar(fill = "#FFFFFF", col = NA)),q2,nrow=2,rel_heights = c(1, 8))
          ,plot_grid(textbox_grob("In-degree centrality", gp = gpar(fontsize = 16, col='black'), background_gp = gpar(fill = "#FFFFFF", col = NA)),q3,nrow=2,rel_heights = c(1, 8)),nrow=3)

cor_rank <- bc_rank[which(bc_rank$cor*2 > bc_rank$all),]%>% arrange(desc(cor))
cor_rank$GeneName <- rownames(cor_rank)

cor_rank

options(repr.plot.width=6, repr.plot.height=4)
ggplot(cor_rank[1:10,], aes(x=reorder(GeneName, cor, decreasing = FALSE),  y=cor)) + geom_point(size=4)+ 
  labs(title="Cortex-specific TF Prioritization",x="", y = "Combined centrality score")+
  theme_classic()+
  theme(axis.text.x = element_text(color = "black", 
                           size = 16, angle = 0),
          axis.text.y = element_text( color = "black", 
                           size = 16, angle = 0),
        plot.title = element_text(color="black", size=18, face="bold"),
        axis.title.x = element_text(color="black", size=16, face="bold"))+ coord_flip()

write.csv(cor_rank,"Cortex_TF_centrality_all_transition_zscore3.csv", quote=FALSE, row.names=TRUE)

tf_rank <- cor_rank %>% rownames(.)
p1 <- plot_bc(tf_rank[1]) + plot_bc(tf_rank[2]) + plot_bc(tf_rank[3]) + plot_bc(tf_rank[4]) + plot_bc(tf_rank[5])
p2 <- plot_bc(tf_rank[6]) + plot_bc(tf_rank[7]) + plot_bc(tf_rank[8]) + plot_bc(tf_rank[9]) + plot_bc(tf_rank[10])
p3 <- plot_bc(tf_rank[11]) + plot_bc(tf_rank[12]) + plot_bc(tf_rank[13]) + plot_bc(tf_rank[14]) + plot_bc(tf_rank[15])
#p4 <- plot_bc(tf_rank[16]) + plot_blank("Empty")+ plot_blank("Empty") + plot_blank("Empty")+ plot_blank("Empty")

p5 <- plot_oc(tf_rank[1]) + plot_oc(tf_rank[2]) + plot_oc(tf_rank[3]) + plot_oc(tf_rank[4]) + plot_oc(tf_rank[5])
p6 <- plot_oc(tf_rank[6]) + plot_oc(tf_rank[7]) + plot_oc(tf_rank[8]) + plot_oc(tf_rank[9]) + plot_oc(tf_rank[10])
p7 <- plot_oc(tf_rank[11]) + plot_oc(tf_rank[12]) + plot_oc(tf_rank[13]) + plot_oc(tf_rank[14]) + plot_oc(tf_rank[15])
#p8 <- plot_oc(tf_rank[16]) + plot_blank("Empty")+ plot_blank("Empty") + plot_blank("Empty")+ plot_blank("Empty")

p9 <- plot_ic(tf_rank[1]) + plot_ic(tf_rank[2]) + plot_ic(tf_rank[3]) + plot_ic(tf_rank[4]) + plot_ic(tf_rank[5])
p10 <- plot_ic(tf_rank[6]) + plot_ic(tf_rank[7]) + plot_ic(tf_rank[8]) + plot_ic(tf_rank[9]) + plot_ic(tf_rank[10])
p11 <- plot_ic(tf_rank[11]) + plot_ic(tf_rank[12]) + plot_ic(tf_rank[13]) + plot_ic(tf_rank[14]) + plot_ic(tf_rank[15])
#p12 <- plot_ic(tf_rank[16]) + plot_blank("Empty")+ plot_blank("Empty") + plot_blank("Empty")+ plot_blank("Empty")

p1 <- grid.grabExpr(draw(p1))
p2 <- grid.grabExpr(draw(p2))
p3 <- grid.grabExpr(draw(p3))
#p4 <- grid.grabExpr(draw(p4))
p5 <- grid.grabExpr(draw(p5))
p6 <- grid.grabExpr(draw(p6))
p7 <- grid.grabExpr(draw(p7))
#p8 <- grid.grabExpr(draw(p8))
p9 <- grid.grabExpr(draw(p9))
p10 <- grid.grabExpr(draw(p10))
p11 <- grid.grabExpr(draw(p11))
#p12 <- grid.grabExpr(draw(p12))
col_fun = colorRamp2(c(0, 0.001, 1), c('navy',"khaki1", "red"))
legend_grob1 <- grid.grabExpr(draw(Legend(col_fun = col_fun, title = "weighted\nnetwork\ncentrality")))
col_fun = colorRamp2(c(0, 0.001, 1), c('gainsboro', "grey", "black"))
legend_grob2 <- grid.grabExpr(draw(Legend(col_fun = col_fun, title = "unweighted\nnetwork\ncentrality")))

options(repr.plot.width=24, repr.plot.height=9)
plot_grid(plot_grid(textbox_grob("Main Regulatory TFs for Cortex ", gp = gpar(fontsize = 20, col='black'), background_gp = gpar(fill = "#FFFFFF", col = NA)),
          plot_grid(plot_grid(textbox_grob("Betweenness centrality", gp = gpar(fontsize = 16, col='black'), background_gp = gpar(fill = "#FFFFFF", col = NA)),plot_grid(p1,p2,p3,nrow=3),nrow=2,rel_heights = c(1, 14))
                    ,plot_grid(textbox_grob("Out-degree centrality", gp = gpar(fontsize = 16, col='black'), background_gp = gpar(fill = "#FFFFFF", col = NA)),plot_grid(p5,p6,p7,nrow=3),nrow=2,rel_heights=c(1,14))
                    ,plot_grid(textbox_grob("In-degree centrality", gp = gpar(fontsize = 16, col='black'), background_gp = gpar(fill = "#FFFFFF", col = NA)),plot_grid(p9,p10,p11,nrow=3),nrow=2,rel_heights=c(1,14)),nrow=1),nrow=2,rel_heights = c(1, 14)),plot_grid(legend_grob1,legend_grob2,nrow=2), rel_widths = c(19, 1))

## Top 10 only
tf_rank <- cor_rank %>% rownames(.)
p1 <- plot_bc(tf_rank[1]) + plot_bc(tf_rank[2]) + plot_bc(tf_rank[3]) + plot_bc(tf_rank[4]) + plot_bc(tf_rank[5])
p2 <- plot_bc(tf_rank[6]) + plot_bc(tf_rank[7]) + plot_bc(tf_rank[8]) + plot_bc(tf_rank[9]) + plot_bc(tf_rank[10])
#p3 <- plot_bc(tf_rank[11]) + plot_bc(tf_rank[12]) + plot_bc(tf_rank[13]) + plot_bc(tf_rank[14]) + plot_bc(tf_rank[15])
#p4 <- plot_bc(tf_rank[16]) + plot_blank("Empty")+ plot_blank("Empty") + plot_blank("Empty")+ plot_blank("Empty")

p5 <- plot_oc(tf_rank[1]) + plot_oc(tf_rank[2]) + plot_oc(tf_rank[3]) + plot_oc(tf_rank[4]) + plot_oc(tf_rank[5])
p6 <- plot_oc(tf_rank[6]) + plot_oc(tf_rank[7]) + plot_oc(tf_rank[8]) + plot_oc(tf_rank[9]) + plot_oc(tf_rank[10])
#p7 <- plot_oc(tf_rank[11]) + plot_oc(tf_rank[12]) + plot_oc(tf_rank[13]) + plot_oc(tf_rank[14]) + plot_oc(tf_rank[15])
#p8 <- plot_oc(tf_rank[16]) + plot_blank("Empty")+ plot_blank("Empty") + plot_blank("Empty")+ plot_blank("Empty")

p9 <- plot_ic(tf_rank[1]) + plot_ic(tf_rank[2]) + plot_ic(tf_rank[3]) + plot_ic(tf_rank[4]) + plot_ic(tf_rank[5])
p10 <- plot_ic(tf_rank[6]) + plot_ic(tf_rank[7]) + plot_ic(tf_rank[8]) + plot_ic(tf_rank[9]) + plot_ic(tf_rank[10])
#p11 <- plot_ic(tf_rank[11]) + plot_ic(tf_rank[12]) + plot_ic(tf_rank[13]) + plot_ic(tf_rank[14]) + plot_ic(tf_rank[15])
#p12 <- plot_ic(tf_rank[16]) + plot_blank("Empty")+ plot_blank("Empty") + plot_blank("Empty")+ plot_blank("Empty")

p1 <- grid.grabExpr(draw(p1))
p2 <- grid.grabExpr(draw(p2))
#p3 <- grid.grabExpr(draw(p3))
#p4 <- grid.grabExpr(draw(p4))
p5 <- grid.grabExpr(draw(p5))
p6 <- grid.grabExpr(draw(p6))
#p7 <- grid.grabExpr(draw(p7))
#p8 <- grid.grabExpr(draw(p8))
p9 <- grid.grabExpr(draw(p9))
p10 <- grid.grabExpr(draw(p10))
#p11 <- grid.grabExpr(draw(p11))
#p12 <- grid.grabExpr(draw(p12))
col_fun = colorRamp2(c(0, 0.001, 1), c('navy',"khaki1", "red"))
legend_grob1 <- grid.grabExpr(draw(Legend(col_fun = col_fun, title = "weighted\nnetwork\ncentrality")))
col_fun = colorRamp2(c(0, 0.001, 1), c('gainsboro', "grey", "black"))
legend_grob2 <- grid.grabExpr(draw(Legend(col_fun = col_fun, title = "unweighted\nnetwork\ncentrality")))

options(repr.plot.width=24, repr.plot.height=5)
plot_grid(plot_grid(textbox_grob("Main Regulatory TFs for Cortex ", gp = gpar(fontsize = 20, col='black'), background_gp = gpar(fill = "#FFFFFF", col = NA)),
          plot_grid(plot_grid(textbox_grob("Betweenness centrality", gp = gpar(fontsize = 16, col='black'), background_gp = gpar(fill = "#FFFFFF", col = NA)),plot_grid(p1,p2,nrow=2),nrow=2,rel_heights = c(1, 14))
                    ,plot_grid(textbox_grob("Out-degree centrality", gp = gpar(fontsize = 16, col='black'), background_gp = gpar(fill = "#FFFFFF", col = NA)),plot_grid(p5,p6,nrow=2),nrow=2,rel_heights=c(1,14))
                    ,plot_grid(textbox_grob("In-degree centrality", gp = gpar(fontsize = 16, col='black'), background_gp = gpar(fill = "#FFFFFF", col = NA)),plot_grid(p9,p10,nrow=2),nrow=2,rel_heights=c(1,14)),nrow=1),nrow=2,rel_heights = c(1, 14)),plot_grid(legend_grob1,legend_grob2,nrow=2), rel_widths = c(19, 1))

q1 <- plot_bc(tf_rank[1]) + plot_bc(tf_rank[2]) + plot_bc(tf_rank[3])
q2 <- plot_oc(tf_rank[1]) + plot_oc(tf_rank[2]) + plot_oc(tf_rank[3])
q3 <- plot_ic(tf_rank[1]) + plot_ic(tf_rank[2]) + plot_ic(tf_rank[3])

q1 <- grid.grabExpr(draw(q1))
q2 <- grid.grabExpr(draw(q2))
q3 <- grid.grabExpr(draw(q3))

options(repr.plot.width=6, repr.plot.height=8)
plot_grid(plot_grid(textbox_grob("Betweenness centrality", gp = gpar(fontsize = 16, col='black'), background_gp = gpar(fill = "#FFFFFF", col = NA)),q1,nrow=2,rel_heights = c(1, 8))
          ,plot_grid(textbox_grob("Out-degree centrality", gp = gpar(fontsize = 16, col='black'), background_gp = gpar(fill = "#FFFFFF", col = NA)),q2,nrow=2,rel_heights = c(1, 8))
          ,plot_grid(textbox_grob("In-degree centrality", gp = gpar(fontsize = 16, col='black'), background_gp = gpar(fill = "#FFFFFF", col = NA)),q3,nrow=2,rel_heights = c(1, 8)),nrow=3)

end_rank <- bc_rank[which(bc_rank$end*2 > bc_rank$all),]%>% arrange(desc(end))
end_rank$GeneName <- rownames(end_rank)

end_rank

options(repr.plot.width=6, repr.plot.height=4)
ggplot(end_rank[1:10,], aes(x=reorder(GeneName, end, decreasing = FALSE),  y=end)) + geom_point(size=4)+ 
  labs(title="Endodermis-specific TF Prioritization",x="", y = "Combined centrality score")+
  theme_classic()+
  theme(axis.text.x = element_text(color = "black", 
                           size = 16, angle = 0),
          axis.text.y = element_text( color = "black", 
                           size = 16, angle = 0),
        plot.title = element_text(color="black", size=18, face="bold"),
        axis.title.x = element_text(color="black", size=16, face="bold"))+ coord_flip()

write.csv(end_rank,"Endodermis_TF_centrality_all_transition_zscore3.csv", quote=FALSE, row.names=TRUE)

tf_rank <- end_rank %>% rownames(.)
# Max 20
p1 <- plot_bc(tf_rank[1]) + plot_bc(tf_rank[2]) + plot_bc(tf_rank[3]) + plot_bc(tf_rank[4]) + plot_bc(tf_rank[5])
p2 <- plot_bc(tf_rank[6]) + plot_bc(tf_rank[7]) + plot_bc(tf_rank[8]) + plot_bc(tf_rank[9]) + plot_bc(tf_rank[10])
p3 <- plot_bc(tf_rank[11]) + plot_bc(tf_rank[12]) + plot_bc(tf_rank[13]) + plot_bc(tf_rank[14]) + plot_bc(tf_rank[15])
p4 <- plot_bc(tf_rank[16]) + plot_bc(tf_rank[17]) + plot_bc(tf_rank[18]) + plot_bc(tf_rank[19]) + plot_bc(tf_rank[20])

p5 <- plot_oc(tf_rank[1]) + plot_oc(tf_rank[2]) + plot_oc(tf_rank[3]) + plot_oc(tf_rank[4]) + plot_oc(tf_rank[5])
p6 <- plot_oc(tf_rank[6]) + plot_oc(tf_rank[7]) + plot_oc(tf_rank[8]) + plot_oc(tf_rank[9]) + plot_oc(tf_rank[10])
p7 <- plot_oc(tf_rank[11]) + plot_oc(tf_rank[12]) + plot_oc(tf_rank[13]) + plot_oc(tf_rank[14]) + plot_oc(tf_rank[15])
p8 <- plot_oc(tf_rank[16]) + plot_oc(tf_rank[17]) + plot_oc(tf_rank[18]) + plot_oc(tf_rank[19]) + plot_oc(tf_rank[20])

p9 <- plot_ic(tf_rank[1]) + plot_ic(tf_rank[2]) + plot_ic(tf_rank[3]) + plot_ic(tf_rank[4]) + plot_ic(tf_rank[5])
p10 <- plot_ic(tf_rank[6]) + plot_ic(tf_rank[7]) + plot_ic(tf_rank[8]) + plot_ic(tf_rank[9]) + plot_ic(tf_rank[10])
p11 <- plot_ic(tf_rank[11]) + plot_ic(tf_rank[12]) + plot_ic(tf_rank[13]) + plot_ic(tf_rank[14]) + plot_ic(tf_rank[15])
p12 <- plot_ic(tf_rank[16]) + plot_ic(tf_rank[17]) + plot_ic(tf_rank[18]) + plot_ic(tf_rank[19]) + plot_ic(tf_rank[20])

p1 <- grid.grabExpr(draw(p1))
p2 <- grid.grabExpr(draw(p2))
p3 <- grid.grabExpr(draw(p3))
p4 <- grid.grabExpr(draw(p4))
p5 <- grid.grabExpr(draw(p5))
p6 <- grid.grabExpr(draw(p6))
p7 <- grid.grabExpr(draw(p7))
p8 <- grid.grabExpr(draw(p8))
p9 <- grid.grabExpr(draw(p9))
p10 <- grid.grabExpr(draw(p10))
p11 <- grid.grabExpr(draw(p11))
p12 <- grid.grabExpr(draw(p12))
col_fun = colorRamp2(c(0, 0.001, 1), c('navy',"khaki1", "red"))
legend_grob1 <- grid.grabExpr(draw(Legend(col_fun = col_fun, title = "weighted\nnetwork\ncentrality")))
col_fun = colorRamp2(c(0, 0.001, 1), c('gainsboro', "grey", "black"))
legend_grob2 <- grid.grabExpr(draw(Legend(col_fun = col_fun, title = "unweighted\nnetwork\ncentrality")))

options(repr.plot.width=24, repr.plot.height=12)
plot_grid(plot_grid(textbox_grob("Main Regulatory TFs for Endodermis ", gp = gpar(fontsize = 20, col='black'), background_gp = gpar(fill = "#FFFFFF", col = NA)),
          plot_grid(plot_grid(textbox_grob("Betweenness centrality", gp = gpar(fontsize = 16, col='black'), background_gp = gpar(fill = "#FFFFFF", col = NA)),plot_grid(p1,p2,p3,p4,nrow=4),nrow=2,rel_heights = c(1, 14))
                    ,plot_grid(textbox_grob("Out-degree centrality", gp = gpar(fontsize = 16, col='black'), background_gp = gpar(fill = "#FFFFFF", col = NA)),plot_grid(p5,p6,p7,p8,nrow=4),nrow=2,rel_heights=c(1,14))
                    ,plot_grid(textbox_grob("In-degree centrality", gp = gpar(fontsize = 16, col='black'), background_gp = gpar(fill = "#FFFFFF", col = NA)),plot_grid(p9,p10,p11,p12,nrow=4),nrow=2,rel_heights=c(1,14)),nrow=1),nrow=2,rel_heights = c(1, 14)),plot_grid(legend_grob1,legend_grob2,nrow=2), rel_widths = c(19, 1))

tf_rank <- end_rank %>% rownames(.)
# Max 10
p1 <- plot_bc(tf_rank[1]) + plot_bc(tf_rank[2]) + plot_bc(tf_rank[3]) + plot_bc(tf_rank[4]) + plot_bc(tf_rank[5])
p2 <- plot_bc(tf_rank[6]) + plot_bc(tf_rank[7]) + plot_bc(tf_rank[8]) + plot_bc(tf_rank[9]) + plot_bc(tf_rank[10])
#p3 <- plot_bc(tf_rank[11]) + plot_bc(tf_rank[12]) + plot_bc(tf_rank[13]) + plot_bc(tf_rank[14]) + plot_bc(tf_rank[15])
#p4 <- plot_bc(tf_rank[16]) + plot_bc(tf_rank[17]) + plot_bc(tf_rank[18]) + plot_bc(tf_rank[19]) + plot_bc(tf_rank[20])

p5 <- plot_oc(tf_rank[1]) + plot_oc(tf_rank[2]) + plot_oc(tf_rank[3]) + plot_oc(tf_rank[4]) + plot_oc(tf_rank[5])
p6 <- plot_oc(tf_rank[6]) + plot_oc(tf_rank[7]) + plot_oc(tf_rank[8]) + plot_oc(tf_rank[9]) + plot_oc(tf_rank[10])
#p7 <- plot_oc(tf_rank[11]) + plot_oc(tf_rank[12]) + plot_oc(tf_rank[13]) + plot_oc(tf_rank[14]) + plot_oc(tf_rank[15])
#p8 <- plot_oc(tf_rank[16]) + plot_oc(tf_rank[17]) + plot_oc(tf_rank[18]) + plot_oc(tf_rank[19]) + plot_oc(tf_rank[20])

p9 <- plot_ic(tf_rank[1]) + plot_ic(tf_rank[2]) + plot_ic(tf_rank[3]) + plot_ic(tf_rank[4]) + plot_ic(tf_rank[5])
p10 <- plot_ic(tf_rank[6]) + plot_ic(tf_rank[7]) + plot_ic(tf_rank[8]) + plot_ic(tf_rank[9]) + plot_ic(tf_rank[10])
#p11 <- plot_ic(tf_rank[11]) + plot_ic(tf_rank[12]) + plot_ic(tf_rank[13]) + plot_ic(tf_rank[14]) + plot_ic(tf_rank[15])
#p12 <- plot_ic(tf_rank[16]) + plot_ic(tf_rank[17]) + plot_ic(tf_rank[18]) + plot_ic(tf_rank[19]) + plot_ic(tf_rank[20])

p1 <- grid.grabExpr(draw(p1))
p2 <- grid.grabExpr(draw(p2))
#p3 <- grid.grabExpr(draw(p3))
#p4 <- grid.grabExpr(draw(p4))
p5 <- grid.grabExpr(draw(p5))
p6 <- grid.grabExpr(draw(p6))
#p7 <- grid.grabExpr(draw(p7))
#p8 <- grid.grabExpr(draw(p8))
p9 <- grid.grabExpr(draw(p9))
p10 <- grid.grabExpr(draw(p10))
#p11 <- grid.grabExpr(draw(p11))
#p12 <- grid.grabExpr(draw(p12))
col_fun = colorRamp2(c(0, 0.001, 1), c('navy',"khaki1", "red"))
legend_grob1 <- grid.grabExpr(draw(Legend(col_fun = col_fun, title = "weighted\nnetwork\ncentrality")))
col_fun = colorRamp2(c(0, 0.001, 1), c('gainsboro', "grey", "black"))
legend_grob2 <- grid.grabExpr(draw(Legend(col_fun = col_fun, title = "unweighted\nnetwork\ncentrality")))

options(repr.plot.width=24, repr.plot.height=5)
plot_grid(plot_grid(textbox_grob("Main Regulatory TFs for Endodermis ", gp = gpar(fontsize = 20, col='black'), background_gp = gpar(fill = "#FFFFFF", col = NA)),
          plot_grid(plot_grid(textbox_grob("Betweenness centrality", gp = gpar(fontsize = 16, col='black'), background_gp = gpar(fill = "#FFFFFF", col = NA)),plot_grid(p1,p2,nrow=2),nrow=2,rel_heights = c(1, 14))
                    ,plot_grid(textbox_grob("Out-degree centrality", gp = gpar(fontsize = 16, col='black'), background_gp = gpar(fill = "#FFFFFF", col = NA)),plot_grid(p5,p6,nrow=2),nrow=2,rel_heights=c(1,14))
                    ,plot_grid(textbox_grob("In-degree centrality", gp = gpar(fontsize = 16, col='black'), background_gp = gpar(fill = "#FFFFFF", col = NA)),plot_grid(p9,p10,nrow=2),nrow=2,rel_heights=c(1,14)),nrow=1),nrow=2,rel_heights = c(1, 14)),plot_grid(legend_grob1,legend_grob2,nrow=2), rel_widths = c(19, 1))

q1 <- plot_bc(tf_rank[1]) + plot_bc(tf_rank[2]) + plot_bc(tf_rank[3]) + plot_bc(tf_rank[4]) + plot_bc(tf_rank[5]) + plot_bc(tf_rank[6]) + plot_bc(tf_rank[7]) + plot_bc(tf_rank[8]) + plot_bc(tf_rank[9]) + plot_bc(tf_rank[10]) 
q2 <- plot_oc(tf_rank[1]) + plot_oc(tf_rank[2]) + plot_oc(tf_rank[3]) + plot_oc(tf_rank[4]) + plot_oc(tf_rank[5]) + plot_oc(tf_rank[6]) + plot_oc(tf_rank[7]) + plot_oc(tf_rank[8]) + plot_oc(tf_rank[9]) + plot_oc(tf_rank[10]) 
q3 <- plot_ic(tf_rank[1]) + plot_ic(tf_rank[2]) + plot_ic(tf_rank[3]) + plot_ic(tf_rank[4]) + plot_ic(tf_rank[5]) + plot_ic(tf_rank[6]) + plot_ic(tf_rank[7]) + plot_ic(tf_rank[8]) + plot_ic(tf_rank[9]) + plot_ic(tf_rank[10]) 

q1 <- grid.grabExpr(draw(q1))
q2 <- grid.grabExpr(draw(q2))
q3 <- grid.grabExpr(draw(q3))

options(repr.plot.width=16, repr.plot.height=8)
plot_grid(plot_grid(textbox_grob("Betweenness centrality", gp = gpar(fontsize = 16, col='black'), background_gp = gpar(fill = "#FFFFFF", col = NA)),q1,nrow=2,rel_heights = c(1, 8))
          ,plot_grid(textbox_grob("Out-degree centrality", gp = gpar(fontsize = 16, col='black'), background_gp = gpar(fill = "#FFFFFF", col = NA)),q2,nrow=2,rel_heights = c(1, 8))
          ,plot_grid(textbox_grob("In-degree centrality", gp = gpar(fontsize = 16, col='black'), background_gp = gpar(fill = "#FFFFFF", col = NA)),q3,nrow=2,rel_heights = c(1, 8)),nrow=3)

per_rank <- bc_rank[which(bc_rank$per*2 > bc_rank$all),]%>% arrange(desc(per))
per_rank$GeneName <- rownames(per_rank)

per_rank

options(repr.plot.width=6, repr.plot.height=4)
ggplot(per_rank[1:10,], aes(x=reorder(GeneName, per, decreasing = FALSE),  y=per)) + geom_point(size=4)+ 
  labs(title="Pericycle-specific TF Prioritization",x="", y = "Combined centrality score")+
  theme_classic()+
  theme(axis.text.x = element_text(color = "black", 
                           size = 16, angle = 0),
          axis.text.y = element_text( color = "black", 
                           size = 16, angle = 0),
        plot.title = element_text(color="black", size=18, face="bold"),
        axis.title.x = element_text(color="black", size=16, face="bold"))+ coord_flip()

write.csv(per_rank,"Pericycle_TF_centrality_all_transition_zscore3.csv", quote=FALSE, row.names=TRUE)

tf_rank <- per_rank %>% rownames(.)
# Max 20
p1 <- plot_bc(tf_rank[1]) + plot_bc(tf_rank[2]) + plot_bc(tf_rank[3]) + plot_bc(tf_rank[4]) + plot_bc(tf_rank[5])
p2 <- plot_bc(tf_rank[6]) + plot_bc(tf_rank[7]) + plot_bc(tf_rank[8]) + plot_bc(tf_rank[9]) + plot_bc(tf_rank[10])
p3 <- plot_bc(tf_rank[11]) + plot_bc(tf_rank[12]) + plot_bc(tf_rank[13]) + plot_bc(tf_rank[14]) + plot_bc(tf_rank[15])
p4 <- plot_bc(tf_rank[16]) + plot_bc(tf_rank[17]) + plot_bc(tf_rank[18]) + plot_bc(tf_rank[19]) + plot_bc(tf_rank[20])

p5 <- plot_oc(tf_rank[1]) + plot_oc(tf_rank[2]) + plot_oc(tf_rank[3]) + plot_oc(tf_rank[4]) + plot_oc(tf_rank[5])
p6 <- plot_oc(tf_rank[6]) + plot_oc(tf_rank[7]) + plot_oc(tf_rank[8]) + plot_oc(tf_rank[9]) + plot_oc(tf_rank[10])
p7 <- plot_oc(tf_rank[11]) + plot_oc(tf_rank[12]) + plot_oc(tf_rank[13]) + plot_oc(tf_rank[14]) + plot_oc(tf_rank[15])
p8 <- plot_oc(tf_rank[16]) + plot_oc(tf_rank[17]) + plot_oc(tf_rank[18]) + plot_oc(tf_rank[19]) + plot_oc(tf_rank[20])

p9 <- plot_ic(tf_rank[1]) + plot_ic(tf_rank[2]) + plot_ic(tf_rank[3]) + plot_ic(tf_rank[4]) + plot_ic(tf_rank[5])
p10 <- plot_ic(tf_rank[6]) + plot_ic(tf_rank[7]) + plot_ic(tf_rank[8]) + plot_ic(tf_rank[9]) + plot_ic(tf_rank[10])
p11 <- plot_ic(tf_rank[11]) + plot_ic(tf_rank[12]) + plot_ic(tf_rank[13]) + plot_ic(tf_rank[14]) + plot_ic(tf_rank[15])
p12 <- plot_ic(tf_rank[16]) + plot_ic(tf_rank[17]) + plot_ic(tf_rank[18]) + plot_ic(tf_rank[19]) + plot_ic(tf_rank[20])

p1 <- grid.grabExpr(draw(p1))
p2 <- grid.grabExpr(draw(p2))
p3 <- grid.grabExpr(draw(p3))
p4 <- grid.grabExpr(draw(p4))
p5 <- grid.grabExpr(draw(p5))
p6 <- grid.grabExpr(draw(p6))
p7 <- grid.grabExpr(draw(p7))
p8 <- grid.grabExpr(draw(p8))
p9 <- grid.grabExpr(draw(p9))
p10 <- grid.grabExpr(draw(p10))
p11 <- grid.grabExpr(draw(p11))
p12 <- grid.grabExpr(draw(p12))
col_fun = colorRamp2(c(0, 0.001, 1), c('navy',"khaki1", "red"))
legend_grob1 <- grid.grabExpr(draw(Legend(col_fun = col_fun, title = "weighted\nnetwork\ncentrality")))
col_fun = colorRamp2(c(0, 0.001, 1), c('gainsboro', "grey", "black"))
legend_grob2 <- grid.grabExpr(draw(Legend(col_fun = col_fun, title = "unweighted\nnetwork\ncentrality")))

options(repr.plot.width=24, repr.plot.height=12)
plot_grid(plot_grid(textbox_grob("Main Regulatory TFs for Pericycle ", gp = gpar(fontsize = 20, col='black'), background_gp = gpar(fill = "#FFFFFF", col = NA)),
          plot_grid(plot_grid(textbox_grob("Betweenness centrality", gp = gpar(fontsize = 16, col='black'), background_gp = gpar(fill = "#FFFFFF", col = NA)),plot_grid(p1,p2,p3,p4,nrow=4),nrow=2,rel_heights = c(1, 14))
                    ,plot_grid(textbox_grob("Out-degree centrality", gp = gpar(fontsize = 16, col='black'), background_gp = gpar(fill = "#FFFFFF", col = NA)),plot_grid(p5,p6,p7,p8,nrow=4),nrow=2,rel_heights=c(1,14))
                    ,plot_grid(textbox_grob("In-degree centrality", gp = gpar(fontsize = 16, col='black'), background_gp = gpar(fill = "#FFFFFF", col = NA)),plot_grid(p9,p10,p11,p12,nrow=4),nrow=2,rel_heights=c(1,14)),nrow=1),nrow=2,rel_heights = c(1, 14)),plot_grid(legend_grob1,legend_grob2,nrow=2), rel_widths = c(19, 1))

q1 <- plot_bc(tf_rank[1]) + plot_bc(tf_rank[2]) + plot_bc(tf_rank[3])
q2 <- plot_oc(tf_rank[1]) + plot_oc(tf_rank[2]) + plot_oc(tf_rank[3])
q3 <- plot_ic(tf_rank[1]) + plot_ic(tf_rank[2]) + plot_ic(tf_rank[3])

q1 <- grid.grabExpr(draw(q1))
q2 <- grid.grabExpr(draw(q2))
q3 <- grid.grabExpr(draw(q3))

options(repr.plot.width=6, repr.plot.height=8)
plot_grid(plot_grid(textbox_grob("Betweenness centrality", gp = gpar(fontsize = 16, col='black'), background_gp = gpar(fill = "#FFFFFF", col = NA)),q1,nrow=2,rel_heights = c(1, 8))
          ,plot_grid(textbox_grob("Out-degree centrality", gp = gpar(fontsize = 16, col='black'), background_gp = gpar(fill = "#FFFFFF", col = NA)),q2,nrow=2,rel_heights = c(1, 8))
          ,plot_grid(textbox_grob("In-degree centrality", gp = gpar(fontsize = 16, col='black'), background_gp = gpar(fill = "#FFFFFF", col = NA)),q3,nrow=2,rel_heights = c(1, 8)),nrow=3)

pro_rank <- bc_rank[which(bc_rank$pro*2 > bc_rank$all),]%>% arrange(desc(pro))
pro_rank$GeneName <- rownames(pro_rank)

pro_rank

options(repr.plot.width=6, repr.plot.height=4)
ggplot(pro_rank[1:10,], aes(x=reorder(GeneName, pro, decreasing = FALSE),  y=pro)) + geom_point(size=4)+ 
  labs(title="Procambium-specific TF Prioritization",x="", y = "Combined centrality score")+
  theme_classic()+
  theme(axis.text.x = element_text(color = "black", 
                           size = 16, angle = 0),
          axis.text.y = element_text( color = "black", 
                           size = 16, angle = 0),
        plot.title = element_text(color="black", size=18, face="bold"),
        axis.title.x = element_text(color="black", size=16, face="bold"))+ coord_flip()

write.csv(pro_rank,"Procambium_TF_centrality_all_transition_zscore3.csv", quote=FALSE, row.names=TRUE)

tf_rank <- pro_rank %>% rownames(.)
# Max 30
p1 <- plot_bc(tf_rank[1]) + plot_bc(tf_rank[2]) + plot_bc(tf_rank[3]) + plot_bc(tf_rank[4]) + plot_bc(tf_rank[5])
p2 <- plot_bc(tf_rank[6]) + plot_bc(tf_rank[7]) + plot_bc(tf_rank[8]) + plot_bc(tf_rank[9]) + plot_bc(tf_rank[10])
p3 <- plot_bc(tf_rank[11]) + plot_bc(tf_rank[12]) + plot_bc(tf_rank[13]) + plot_bc(tf_rank[14]) + plot_blank(tf_rank[15])

p4 <- plot_oc(tf_rank[1]) + plot_oc(tf_rank[2]) + plot_oc(tf_rank[3]) + plot_oc(tf_rank[4]) + plot_oc(tf_rank[5])
p5 <- plot_oc(tf_rank[6]) + plot_oc(tf_rank[7]) + plot_oc(tf_rank[8]) + plot_oc(tf_rank[9]) + plot_oc(tf_rank[10])
p6 <- plot_oc(tf_rank[11]) + plot_oc(tf_rank[12]) + plot_oc(tf_rank[13]) + plot_oc(tf_rank[14]) + plot_blank(tf_rank[15])

p7 <- plot_ic(tf_rank[1]) + plot_ic(tf_rank[2]) + plot_ic(tf_rank[3]) + plot_ic(tf_rank[4]) + plot_ic(tf_rank[5])
p8 <- plot_ic(tf_rank[6]) + plot_ic(tf_rank[7]) + plot_ic(tf_rank[8]) + plot_ic(tf_rank[9]) + plot_ic(tf_rank[10])
p9 <- plot_ic(tf_rank[11]) + plot_ic(tf_rank[12]) + plot_ic(tf_rank[13]) + plot_ic(tf_rank[14]) + plot_blank(tf_rank[15])

p1 <- grid.grabExpr(draw(p1))
p2 <- grid.grabExpr(draw(p2))
p3 <- grid.grabExpr(draw(p3))
p4 <- grid.grabExpr(draw(p4))
p5 <- grid.grabExpr(draw(p5))
p6 <- grid.grabExpr(draw(p6))
p7 <- grid.grabExpr(draw(p7))
p8 <- grid.grabExpr(draw(p8))
p9 <- grid.grabExpr(draw(p9))
col_fun = colorRamp2(c(0, 0.001, 1), c('navy',"khaki1", "red"))
legend_grob1 <- grid.grabExpr(draw(Legend(col_fun = col_fun, title = "weighted\nnetwork\ncentrality")))
col_fun = colorRamp2(c(0, 0.001, 1), c('gainsboro', "grey", "black"))
legend_grob2 <- grid.grabExpr(draw(Legend(col_fun = col_fun, title = "unweighted\nnetwork\ncentrality")))

options(repr.plot.width=24, repr.plot.height=10)
plot_grid(plot_grid(textbox_grob("Main Regulatory TFs for Procambium ", gp = gpar(fontsize = 20, col='black'), background_gp = gpar(fill = "#FFFFFF", col = NA)),
          plot_grid(plot_grid(textbox_grob("Betweenness centrality", gp = gpar(fontsize = 16, col='black'), background_gp = gpar(fill = "#FFFFFF", col = NA)),plot_grid(p1,p2,p3,nrow=3),nrow=2,rel_heights = c(1, 14))
                    ,plot_grid(textbox_grob("Out-degree centrality", gp = gpar(fontsize = 16, col='black'), background_gp = gpar(fill = "#FFFFFF", col = NA)),plot_grid(p4,p5,p6,nrow=3),nrow=2,rel_heights=c(1,14))
                    ,plot_grid(textbox_grob("In-degree centrality", gp = gpar(fontsize = 16, col='black'), background_gp = gpar(fill = "#FFFFFF", col = NA)),plot_grid(p7,p8,p9,nrow=3),nrow=2,rel_heights=c(1,14)),nrow=1),nrow=2,rel_heights = c(1, 14)),plot_grid(legend_grob1,legend_grob2,nrow=2), rel_widths = c(19, 1))

phl_rank <- bc_rank[which(bc_rank$phl*2 > bc_rank$all),]%>% arrange(desc(phl))
phl_rank$GeneName <- rownames(phl_rank)

phl_rank

options(repr.plot.width=6, repr.plot.height=4)
ggplot(phl_rank[1:10,], aes(x=reorder(GeneName, phl, decreasing = FALSE),  y=phl)) + geom_point(size=4)+ 
  labs(title="Phloem-specific TF Prioritization",x="", y = "Combined centrality score")+
  theme_classic()+
  theme(axis.text.x = element_text(color = "black", 
                           size = 16, angle = 0),
          axis.text.y = element_text( color = "black", 
                           size = 16, angle = 0),
        plot.title = element_text(color="black", size=18, face="bold"),
        axis.title.x = element_text(color="black", size=16, face="bold"))+ coord_flip()

write.csv(phl_rank,"Phloem_TF_centrality_all_transition_zscore3.csv", quote=FALSE, row.names=TRUE)

tf_rank <- phl_rank %>% rownames(.)
# Max 20
p1 <- plot_bc(tf_rank[1]) + plot_bc(tf_rank[2]) + plot_bc(tf_rank[3]) + plot_bc(tf_rank[4]) + plot_bc(tf_rank[5])
p2 <- plot_bc(tf_rank[6]) + plot_bc(tf_rank[7]) + plot_bc(tf_rank[8]) + plot_bc(tf_rank[9]) + plot_bc(tf_rank[10])
p3 <- plot_bc(tf_rank[11]) + plot_bc(tf_rank[12]) + plot_bc(tf_rank[13]) + plot_bc(tf_rank[14]) + plot_bc(tf_rank[15])
p4 <- plot_bc(tf_rank[16]) + plot_bc(tf_rank[17]) + plot_bc(tf_rank[18]) + plot_bc(tf_rank[19]) + plot_bc(tf_rank[20])

p5 <- plot_oc(tf_rank[1]) + plot_oc(tf_rank[2]) + plot_oc(tf_rank[3]) + plot_oc(tf_rank[4]) + plot_oc(tf_rank[5])
p6 <- plot_oc(tf_rank[6]) + plot_oc(tf_rank[7]) + plot_oc(tf_rank[8]) + plot_oc(tf_rank[9]) + plot_oc(tf_rank[10])
p7 <- plot_oc(tf_rank[11]) + plot_oc(tf_rank[12]) + plot_oc(tf_rank[13]) + plot_oc(tf_rank[14]) + plot_oc(tf_rank[15])
p8 <- plot_oc(tf_rank[16]) + plot_oc(tf_rank[17]) + plot_oc(tf_rank[18]) + plot_oc(tf_rank[19]) + plot_oc(tf_rank[20])

p9 <- plot_ic(tf_rank[1]) + plot_ic(tf_rank[2]) + plot_ic(tf_rank[3]) + plot_ic(tf_rank[4]) + plot_ic(tf_rank[5])
p10 <- plot_ic(tf_rank[6]) + plot_ic(tf_rank[7]) + plot_ic(tf_rank[8]) + plot_ic(tf_rank[9]) + plot_ic(tf_rank[10])
p11 <- plot_ic(tf_rank[11]) + plot_ic(tf_rank[12]) + plot_ic(tf_rank[13]) + plot_ic(tf_rank[14]) + plot_ic(tf_rank[15])
p12 <- plot_ic(tf_rank[16]) + plot_ic(tf_rank[17]) + plot_ic(tf_rank[18]) + plot_ic(tf_rank[19]) + plot_ic(tf_rank[20])

p1 <- grid.grabExpr(draw(p1))
p2 <- grid.grabExpr(draw(p2))
p3 <- grid.grabExpr(draw(p3))
p4 <- grid.grabExpr(draw(p4))
p5 <- grid.grabExpr(draw(p5))
p6 <- grid.grabExpr(draw(p6))
p7 <- grid.grabExpr(draw(p7))
p8 <- grid.grabExpr(draw(p8))
p9 <- grid.grabExpr(draw(p9))
p10 <- grid.grabExpr(draw(p10))
p11 <- grid.grabExpr(draw(p11))
p12 <- grid.grabExpr(draw(p12))
col_fun = colorRamp2(c(0, 0.001, 1), c('navy',"khaki1", "red"))
legend_grob1 <- grid.grabExpr(draw(Legend(col_fun = col_fun, title = "weighted\nnetwork\ncentrality")))
col_fun = colorRamp2(c(0, 0.001, 1), c('gainsboro', "grey", "black"))
legend_grob2 <- grid.grabExpr(draw(Legend(col_fun = col_fun, title = "unweighted\nnetwork\ncentrality")))

options(repr.plot.width=24, repr.plot.height=12)
plot_grid(plot_grid(textbox_grob("Main Regulatory TFs for Phloem ", gp = gpar(fontsize = 20, col='black'), background_gp = gpar(fill = "#FFFFFF", col = NA)),
          plot_grid(plot_grid(textbox_grob("Betweenness centrality", gp = gpar(fontsize = 16, col='black'), background_gp = gpar(fill = "#FFFFFF", col = NA)),plot_grid(p1,p2,p3,p4,nrow=4),nrow=2,rel_heights = c(1, 14))
                    ,plot_grid(textbox_grob("Out-degree centrality", gp = gpar(fontsize = 16, col='black'), background_gp = gpar(fill = "#FFFFFF", col = NA)),plot_grid(p5,p6,p7,p8,nrow=4),nrow=2,rel_heights=c(1,14))
                    ,plot_grid(textbox_grob("In-degree centrality", gp = gpar(fontsize = 16, col='black'), background_gp = gpar(fill = "#FFFFFF", col = NA)),plot_grid(p9,p10,p11,p12,nrow=4),nrow=2,rel_heights=c(1,14)),nrow=1),nrow=2,rel_heights = c(1, 14)),plot_grid(legend_grob1,legend_grob2,nrow=2), rel_widths = c(19, 1))

xyl_rank <- bc_rank[which(bc_rank$xyl*2 > bc_rank$all),]%>% arrange(desc(xyl))
xyl_rank$GeneName <- rownames(xyl_rank)

xyl_rank

options(repr.plot.width=6, repr.plot.height=4)
ggplot(xyl_rank[1:10,], aes(x=reorder(GeneName, xyl, decreasing = FALSE),  y=xyl)) + geom_point(size=4)+ 
  labs(title="Xylem-specific TF Prioritization",x="", y = "Combined centrality score")+
  theme_classic()+
  theme(axis.text.x = element_text(color = "black", 
                           size = 16, angle = 0),
          axis.text.y = element_text( color = "black", 
                           size = 16, angle = 0),
        plot.title = element_text(color="black", size=18, face="bold"),
        axis.title.x = element_text(color="black", size=16, face="bold"))+ coord_flip()

write.csv(xyl_rank,"Xylem_TF_centrality_all_transition_zscore3.csv", quote=FALSE, row.names=TRUE)

tf_rank <- xyl_rank %>% rownames(.)
# Max 20
p1 <- plot_bc(tf_rank[1]) + plot_bc(tf_rank[2]) + plot_bc(tf_rank[3]) + plot_bc(tf_rank[4]) + plot_bc(tf_rank[5])
p2 <- plot_bc(tf_rank[6]) + plot_bc(tf_rank[7]) + plot_bc(tf_rank[8]) + plot_bc(tf_rank[9]) + plot_bc(tf_rank[10])
p3 <- plot_bc(tf_rank[11]) + plot_bc(tf_rank[12]) + plot_bc(tf_rank[13]) + plot_bc(tf_rank[14]) + plot_bc(tf_rank[15])
p4 <- plot_bc(tf_rank[16]) + plot_bc(tf_rank[17]) + plot_bc(tf_rank[18]) + plot_bc(tf_rank[19]) + plot_bc(tf_rank[20])

p5 <- plot_oc(tf_rank[1]) + plot_oc(tf_rank[2]) + plot_oc(tf_rank[3]) + plot_oc(tf_rank[4]) + plot_oc(tf_rank[5])
p6 <- plot_oc(tf_rank[6]) + plot_oc(tf_rank[7]) + plot_oc(tf_rank[8]) + plot_oc(tf_rank[9]) + plot_oc(tf_rank[10])
p7 <- plot_oc(tf_rank[11]) + plot_oc(tf_rank[12]) + plot_oc(tf_rank[13]) + plot_oc(tf_rank[14]) + plot_oc(tf_rank[15])
p8 <- plot_oc(tf_rank[16]) + plot_oc(tf_rank[17]) + plot_oc(tf_rank[18]) + plot_oc(tf_rank[19]) + plot_oc(tf_rank[20])

p9 <- plot_ic(tf_rank[1]) + plot_ic(tf_rank[2]) + plot_ic(tf_rank[3]) + plot_ic(tf_rank[4]) + plot_ic(tf_rank[5])
p10 <- plot_ic(tf_rank[6]) + plot_ic(tf_rank[7]) + plot_ic(tf_rank[8]) + plot_ic(tf_rank[9]) + plot_ic(tf_rank[10])
p11 <- plot_ic(tf_rank[11]) + plot_ic(tf_rank[12]) + plot_ic(tf_rank[13]) + plot_ic(tf_rank[14]) + plot_ic(tf_rank[15])
p12 <- plot_ic(tf_rank[16]) + plot_ic(tf_rank[17]) + plot_ic(tf_rank[18]) + plot_ic(tf_rank[19]) + plot_ic(tf_rank[20])

p1 <- grid.grabExpr(draw(p1))
p2 <- grid.grabExpr(draw(p2))
p3 <- grid.grabExpr(draw(p3))
p4 <- grid.grabExpr(draw(p4))
p5 <- grid.grabExpr(draw(p5))
p6 <- grid.grabExpr(draw(p6))
p7 <- grid.grabExpr(draw(p7))
p8 <- grid.grabExpr(draw(p8))
p9 <- grid.grabExpr(draw(p9))
p10 <- grid.grabExpr(draw(p10))
p11 <- grid.grabExpr(draw(p11))
p12 <- grid.grabExpr(draw(p12))
col_fun = colorRamp2(c(0, 0.001, 1), c('navy',"khaki1", "red"))
legend_grob1 <- grid.grabExpr(draw(Legend(col_fun = col_fun, title = "weighted\nnetwork\ncentrality")))
col_fun = colorRamp2(c(0, 0.001, 1), c('gainsboro', "grey", "black"))
legend_grob2 <- grid.grabExpr(draw(Legend(col_fun = col_fun, title = "unweighted\nnetwork\ncentrality")))

options(repr.plot.width=24, repr.plot.height=12)
plot_grid(plot_grid(textbox_grob("Main Regulatory TFs for Xylem ", gp = gpar(fontsize = 20, col='black'), background_gp = gpar(fill = "#FFFFFF", col = NA)),
          plot_grid(plot_grid(textbox_grob("Betweenness centrality", gp = gpar(fontsize = 16, col='black'), background_gp = gpar(fill = "#FFFFFF", col = NA)),plot_grid(p1,p2,p3,p4,nrow=4),nrow=2,rel_heights = c(1, 14))
                    ,plot_grid(textbox_grob("Out-degree centrality", gp = gpar(fontsize = 16, col='black'), background_gp = gpar(fill = "#FFFFFF", col = NA)),plot_grid(p5,p6,p7,p8,nrow=4),nrow=2,rel_heights=c(1,14))
                    ,plot_grid(textbox_grob("In-degree centrality", gp = gpar(fontsize = 16, col='black'), background_gp = gpar(fill = "#FFFFFF", col = NA)),plot_grid(p9,p10,p11,p12,nrow=4),nrow=2,rel_heights=c(1,14)),nrow=1),nrow=2,rel_heights = c(1, 14)),plot_grid(legend_grob1,legend_grob2,nrow=2), rel_widths = c(19, 1))

lrc_rank <- bc_rank[which(bc_rank$lrc*2 > bc_rank$all),]%>% arrange(desc(lrc))
lrc_rank$GeneName <- rownames(lrc_rank)

lrc_rank

options(repr.plot.width=6, repr.plot.height=4)
ggplot(lrc_rank[1:10,], aes(x=reorder(GeneName, lrc, decreasing = FALSE),  y=lrc)) + geom_point(size=4)+ 
  labs(title="LRC-specific TF Prioritization",x="", y = "Combined centrality score")+
  theme_classic()+
  theme(axis.text.x = element_text(color = "black", 
                           size = 16, angle = 0),
          axis.text.y = element_text( color = "black", 
                           size = 16, angle = 0),
        plot.title = element_text(color="black", size=18, face="bold"),
        axis.title.x = element_text(color="black", size=16, face="bold"))+ coord_flip()

write.csv(lrc_rank,"LRC_TF_centrality_all_transition_zscore3.csv", quote=FALSE, row.names=TRUE)

tf_rank <- lrc_rank %>% rownames(.)
# Max 20
p1 <- plot_bc(tf_rank[1]) + plot_bc(tf_rank[2]) + plot_bc(tf_rank[3]) + plot_bc(tf_rank[4]) + plot_bc(tf_rank[5])
p2 <- plot_bc(tf_rank[6]) + plot_bc(tf_rank[7]) + plot_bc(tf_rank[8]) + plot_bc(tf_rank[9]) + plot_bc(tf_rank[10])
p3 <- plot_bc(tf_rank[11]) + plot_bc(tf_rank[12]) + plot_bc(tf_rank[13]) + plot_bc(tf_rank[14]) + plot_bc(tf_rank[15])
p4 <- plot_bc(tf_rank[16]) + plot_bc(tf_rank[17]) + plot_bc(tf_rank[18]) + plot_bc(tf_rank[19]) + plot_bc(tf_rank[20])

p5 <- plot_oc(tf_rank[1]) + plot_oc(tf_rank[2]) + plot_oc(tf_rank[3]) + plot_oc(tf_rank[4]) + plot_oc(tf_rank[5])
p6 <- plot_oc(tf_rank[6]) + plot_oc(tf_rank[7]) + plot_oc(tf_rank[8]) + plot_oc(tf_rank[9]) + plot_oc(tf_rank[10])
p7 <- plot_oc(tf_rank[11]) + plot_oc(tf_rank[12]) + plot_oc(tf_rank[13]) + plot_oc(tf_rank[14]) + plot_oc(tf_rank[15])
p8 <- plot_oc(tf_rank[16]) + plot_oc(tf_rank[17]) + plot_oc(tf_rank[18]) + plot_oc(tf_rank[19]) + plot_oc(tf_rank[20])

p9 <- plot_ic(tf_rank[1]) + plot_ic(tf_rank[2]) + plot_ic(tf_rank[3]) + plot_ic(tf_rank[4]) + plot_ic(tf_rank[5])
p10 <- plot_ic(tf_rank[6]) + plot_ic(tf_rank[7]) + plot_ic(tf_rank[8]) + plot_ic(tf_rank[9]) + plot_ic(tf_rank[10])
p11 <- plot_ic(tf_rank[11]) + plot_ic(tf_rank[12]) + plot_ic(tf_rank[13]) + plot_ic(tf_rank[14]) + plot_ic(tf_rank[15])
p12 <- plot_ic(tf_rank[16]) + plot_ic(tf_rank[17]) + plot_ic(tf_rank[18]) + plot_ic(tf_rank[19]) + plot_ic(tf_rank[20])

p1 <- grid.grabExpr(draw(p1))
p2 <- grid.grabExpr(draw(p2))
p3 <- grid.grabExpr(draw(p3))
p4 <- grid.grabExpr(draw(p4))
p5 <- grid.grabExpr(draw(p5))
p6 <- grid.grabExpr(draw(p6))
p7 <- grid.grabExpr(draw(p7))
p8 <- grid.grabExpr(draw(p8))
p9 <- grid.grabExpr(draw(p9))
p10 <- grid.grabExpr(draw(p10))
p11 <- grid.grabExpr(draw(p11))
p12 <- grid.grabExpr(draw(p12))
col_fun = colorRamp2(c(0, 0.001, 1), c('navy',"khaki1", "red"))
legend_grob1 <- grid.grabExpr(draw(Legend(col_fun = col_fun, title = "weighted\nnetwork\ncentrality")))
col_fun = colorRamp2(c(0, 0.001, 1), c('gainsboro', "grey", "black"))
legend_grob2 <- grid.grabExpr(draw(Legend(col_fun = col_fun, title = "unweighted\nnetwork\ncentrality")))

options(repr.plot.width=24, repr.plot.height=12)
plot_grid(plot_grid(textbox_grob("Main Regulatory TFs for Lateral Root Cap ", gp = gpar(fontsize = 20, col='black'), background_gp = gpar(fill = "#FFFFFF", col = NA)),
          plot_grid(plot_grid(textbox_grob("Betweenness centrality", gp = gpar(fontsize = 16, col='black'), background_gp = gpar(fill = "#FFFFFF", col = NA)),plot_grid(p1,p2,p3,p4,nrow=4),nrow=2,rel_heights = c(1, 14))
                    ,plot_grid(textbox_grob("Out-degree centrality", gp = gpar(fontsize = 16, col='black'), background_gp = gpar(fill = "#FFFFFF", col = NA)),plot_grid(p5,p6,p7,p8,nrow=4),nrow=2,rel_heights=c(1,14))
                    ,plot_grid(textbox_grob("In-degree centrality", gp = gpar(fontsize = 16, col='black'), background_gp = gpar(fill = "#FFFFFF", col = NA)),plot_grid(p9,p10,p11,p12,nrow=4),nrow=2,rel_heights=c(1,14)),nrow=1),nrow=2,rel_heights = c(1, 14)),plot_grid(legend_grob1,legend_grob2,nrow=2), rel_widths = c(19, 1))

col_rank <- bc_rank[which(bc_rank$col*2 > bc_rank$all),]%>% arrange(desc(col))
col_rank$GeneName <- rownames(col_rank)

col_rank

options(repr.plot.width=6, repr.plot.height=4)
ggplot(col_rank[1:10,], aes(x=reorder(GeneName, col, decreasing = FALSE),  y=col)) + geom_point(size=4)+ 
  labs(title="Columella-specific TF Prioritization",x="", y = "Combined centrality score")+
  theme_classic()+
  theme(axis.text.x = element_text(color = "black", 
                           size = 16, angle = 0),
          axis.text.y = element_text( color = "black", 
                           size = 16, angle = 0),
        plot.title = element_text(color="black", size=18, face="bold"),
        axis.title.x = element_text(color="black", size=16, face="bold"))+ coord_flip()

write.csv(col_rank,"Columella_TF_centrality_all_transition_zscore3.csv", quote=FALSE, row.names=TRUE)

tf_rank <- col_rank %>% rownames(.)
# Max 20
p1 <- plot_bc(tf_rank[1]) + plot_bc(tf_rank[2]) + plot_bc(tf_rank[3]) + plot_bc(tf_rank[4]) + plot_bc(tf_rank[5])
p2 <- plot_bc(tf_rank[6]) + plot_bc(tf_rank[7]) + plot_bc(tf_rank[8]) + plot_bc(tf_rank[9]) + plot_bc(tf_rank[10])
p3 <- plot_bc(tf_rank[11]) + plot_bc(tf_rank[12]) + plot_bc(tf_rank[13]) + plot_bc(tf_rank[14]) + plot_bc(tf_rank[15])
p4 <- plot_bc(tf_rank[16]) + plot_bc(tf_rank[17]) + plot_bc(tf_rank[18]) + plot_bc(tf_rank[19]) + plot_bc(tf_rank[20])

p5 <- plot_oc(tf_rank[1]) + plot_oc(tf_rank[2]) + plot_oc(tf_rank[3]) + plot_oc(tf_rank[4]) + plot_oc(tf_rank[5])
p6 <- plot_oc(tf_rank[6]) + plot_oc(tf_rank[7]) + plot_oc(tf_rank[8]) + plot_oc(tf_rank[9]) + plot_oc(tf_rank[10])
p7 <- plot_oc(tf_rank[11]) + plot_oc(tf_rank[12]) + plot_oc(tf_rank[13]) + plot_oc(tf_rank[14]) + plot_oc(tf_rank[15])
p8 <- plot_oc(tf_rank[16]) + plot_oc(tf_rank[17]) + plot_oc(tf_rank[18]) + plot_oc(tf_rank[19]) + plot_oc(tf_rank[20])

p9 <- plot_ic(tf_rank[1]) + plot_ic(tf_rank[2]) + plot_ic(tf_rank[3]) + plot_ic(tf_rank[4]) + plot_ic(tf_rank[5])
p10 <- plot_ic(tf_rank[6]) + plot_ic(tf_rank[7]) + plot_ic(tf_rank[8]) + plot_ic(tf_rank[9]) + plot_ic(tf_rank[10])
p11 <- plot_ic(tf_rank[11]) + plot_ic(tf_rank[12]) + plot_ic(tf_rank[13]) + plot_ic(tf_rank[14]) + plot_ic(tf_rank[15])
p12 <- plot_ic(tf_rank[16]) + plot_ic(tf_rank[17]) + plot_ic(tf_rank[18]) + plot_ic(tf_rank[19]) + plot_ic(tf_rank[20])

p1 <- grid.grabExpr(draw(p1))
p2 <- grid.grabExpr(draw(p2))
p3 <- grid.grabExpr(draw(p3))
p4 <- grid.grabExpr(draw(p4))
p5 <- grid.grabExpr(draw(p5))
p6 <- grid.grabExpr(draw(p6))
p7 <- grid.grabExpr(draw(p7))
p8 <- grid.grabExpr(draw(p8))
p9 <- grid.grabExpr(draw(p9))
p10 <- grid.grabExpr(draw(p10))
p11 <- grid.grabExpr(draw(p11))
p12 <- grid.grabExpr(draw(p12))
col_fun = colorRamp2(c(0, 0.001, 1), c('navy',"khaki1", "red"))
legend_grob1 <- grid.grabExpr(draw(Legend(col_fun = col_fun, title = "weighted\nnetwork\ncentrality")))
col_fun = colorRamp2(c(0, 0.001, 1), c('gainsboro', "grey", "black"))
legend_grob2 <- grid.grabExpr(draw(Legend(col_fun = col_fun, title = "unweighted\nnetwork\ncentrality")))

options(repr.plot.width=24, repr.plot.height=12)
plot_grid(plot_grid(textbox_grob("Main Regulatory TFs for Columella ", gp = gpar(fontsize = 20, col='black'), background_gp = gpar(fill = "#FFFFFF", col = NA)),
          plot_grid(plot_grid(textbox_grob("Betweenness centrality", gp = gpar(fontsize = 16, col='black'), background_gp = gpar(fill = "#FFFFFF", col = NA)),plot_grid(p1,p2,p3,p4,nrow=4),nrow=2,rel_heights = c(1, 14))
                    ,plot_grid(textbox_grob("Out-degree centrality", gp = gpar(fontsize = 16, col='black'), background_gp = gpar(fill = "#FFFFFF", col = NA)),plot_grid(p5,p6,p7,p8,nrow=4),nrow=2,rel_heights=c(1,14))
                    ,plot_grid(textbox_grob("In-degree centrality", gp = gpar(fontsize = 16, col='black'), background_gp = gpar(fill = "#FFFFFF", col = NA)),plot_grid(p9,p10,p11,p12,nrow=4),nrow=2,rel_heights=c(1,14)),nrow=1),nrow=2,rel_heights = c(1, 14)),plot_grid(legend_grob1,legend_grob2,nrow=2), rel_widths = c(19, 1))

bc_rank <- bc_rank %>% mutate(ground=cor+end, epi=atri+tri, stele=per+pro+xyl+phl, epilrc=atri+tri+lrc, rc=lrc+col)

head(bc_rank)

ground_rank <- bc_rank[which(bc_rank$ground*2 > bc_rank$all),]%>% arrange(desc(ground))
ground_rank <- ground_rank[-c(match(rownames(cor_rank), rownames(ground_rank)),match(rownames(end_rank), rownames(ground_rank))),]
ground_rank$GeneName <- rownames(ground_rank)

ground_rank

options(repr.plot.width=6, repr.plot.height=4)
ggplot(ground_rank[1:10,], aes(x=reorder(GeneName, ground, decreasing = FALSE),  y=ground)) + geom_point(size=4)+ 
  labs(title="Ground Tissue-specific TF Prioritization",x="", y = "Combined centrality score")+
  theme_classic()+
  theme(axis.text.x = element_text(color = "black", 
                           size = 16, angle = 0),
          axis.text.y = element_text( color = "black", 
                           size = 16, angle = 0),
        plot.title = element_text(color="black", size=18, face="bold"),
        axis.title.x = element_text(color="black", size=16, face="bold"))+ coord_flip()

write.csv(ground_rank,"Ground_Tissue_TF_centrality_all_transition_zscore3.csv", quote=FALSE, row.names=TRUE)

tf_rank <- ground_rank %>% rownames(.)
# Max 20
p1 <- plot_bc(tf_rank[1]) + plot_bc(tf_rank[2]) + plot_bc(tf_rank[3]) + plot_bc(tf_rank[4]) + plot_bc(tf_rank[5])
p2 <- plot_bc(tf_rank[6]) + plot_bc(tf_rank[7]) + plot_bc(tf_rank[8]) + plot_bc(tf_rank[9]) + plot_bc(tf_rank[10])
p3 <- plot_bc(tf_rank[11]) + plot_bc(tf_rank[12]) + plot_bc(tf_rank[13]) + plot_bc(tf_rank[14]) + plot_bc(tf_rank[15])
p4 <- plot_bc(tf_rank[16]) + plot_bc(tf_rank[17]) + plot_bc(tf_rank[18]) + plot_bc(tf_rank[19]) + plot_bc(tf_rank[20])

p5 <- plot_oc(tf_rank[1]) + plot_oc(tf_rank[2]) + plot_oc(tf_rank[3]) + plot_oc(tf_rank[4]) + plot_oc(tf_rank[5])
p6 <- plot_oc(tf_rank[6]) + plot_oc(tf_rank[7]) + plot_oc(tf_rank[8]) + plot_oc(tf_rank[9]) + plot_oc(tf_rank[10])
p7 <- plot_oc(tf_rank[11]) + plot_oc(tf_rank[12]) + plot_oc(tf_rank[13]) + plot_oc(tf_rank[14]) + plot_oc(tf_rank[15])
p8 <- plot_oc(tf_rank[16]) + plot_oc(tf_rank[17]) + plot_oc(tf_rank[18]) + plot_oc(tf_rank[19]) + plot_oc(tf_rank[20])

p9 <- plot_ic(tf_rank[1]) + plot_ic(tf_rank[2]) + plot_ic(tf_rank[3]) + plot_ic(tf_rank[4]) + plot_ic(tf_rank[5])
p10 <- plot_ic(tf_rank[6]) + plot_ic(tf_rank[7]) + plot_ic(tf_rank[8]) + plot_ic(tf_rank[9]) + plot_ic(tf_rank[10])
p11 <- plot_ic(tf_rank[11]) + plot_ic(tf_rank[12]) + plot_ic(tf_rank[13]) + plot_ic(tf_rank[14]) + plot_ic(tf_rank[15])
p12 <- plot_ic(tf_rank[16]) + plot_ic(tf_rank[17]) + plot_ic(tf_rank[18]) + plot_ic(tf_rank[19]) + plot_ic(tf_rank[20])

p1 <- grid.grabExpr(draw(p1))
p2 <- grid.grabExpr(draw(p2))
p3 <- grid.grabExpr(draw(p3))
p4 <- grid.grabExpr(draw(p4))
p5 <- grid.grabExpr(draw(p5))
p6 <- grid.grabExpr(draw(p6))
p7 <- grid.grabExpr(draw(p7))
p8 <- grid.grabExpr(draw(p8))
p9 <- grid.grabExpr(draw(p9))
p10 <- grid.grabExpr(draw(p10))
p11 <- grid.grabExpr(draw(p11))
p12 <- grid.grabExpr(draw(p12))
col_fun = colorRamp2(c(0, 0.001, 1), c('navy',"khaki1", "red"))
legend_grob1 <- grid.grabExpr(draw(Legend(col_fun = col_fun, title = "weighted\nnetwork\ncentrality")))
col_fun = colorRamp2(c(0, 0.001, 1), c('gainsboro', "grey", "black"))
legend_grob2 <- grid.grabExpr(draw(Legend(col_fun = col_fun, title = "unweighted\nnetwork\ncentrality")))

options(repr.plot.width=24, repr.plot.height=12)
plot_grid(plot_grid(textbox_grob("Main Regulatory TFs for Ground Tissue ", gp = gpar(fontsize = 20, col='black'), background_gp = gpar(fill = "#FFFFFF", col = NA)),
          plot_grid(plot_grid(textbox_grob("Betweenness centrality", gp = gpar(fontsize = 16, col='black'), background_gp = gpar(fill = "#FFFFFF", col = NA)),plot_grid(p1,p2,p3,p4,nrow=4),nrow=2,rel_heights = c(1, 14))
                    ,plot_grid(textbox_grob("Out-degree centrality", gp = gpar(fontsize = 16, col='black'), background_gp = gpar(fill = "#FFFFFF", col = NA)),plot_grid(p5,p6,p7,p8,nrow=4),nrow=2,rel_heights=c(1,14))
                    ,plot_grid(textbox_grob("In-degree centrality", gp = gpar(fontsize = 16, col='black'), background_gp = gpar(fill = "#FFFFFF", col = NA)),plot_grid(p9,p10,p11,p12,nrow=4),nrow=2,rel_heights=c(1,14)),nrow=1),nrow=2,rel_heights = c(1, 14)),plot_grid(legend_grob1,legend_grob2,nrow=2), rel_widths = c(19, 1))

epi_rank <- bc_rank[which(bc_rank$epi*2 > bc_rank$all),]%>% arrange(desc(epi))
epi_rank <- epi_rank[-c(match(rownames(atri_rank), rownames(epi_rank)),match(rownames(tri_rank), rownames(epi_rank))),]
epi_rank$GeneName <- rownames(epi_rank)

epi_rank

options(repr.plot.width=6, repr.plot.height=4)
ggplot(epi_rank[1:10,], aes(x=reorder(GeneName, epi, decreasing = FALSE),  y=epi)) + geom_point(size=4)+ 
  labs(title="Epidermis-specific TF Prioritization",x="", y = "Combined centrality score")+
  theme_classic()+
  theme(axis.text.x = element_text(color = "black", 
                           size = 16, angle = 0),
          axis.text.y = element_text( color = "black", 
                           size = 16, angle = 0),
        plot.title = element_text(color="black", size=18, face="bold"),
        axis.title.x = element_text(color="black", size=16, face="bold"))+ coord_flip()

write.csv(epi_rank,"Epidermis_TF_centrality_all_transition_zscore3.csv", quote=FALSE, row.names=TRUE)

tf_rank <- epi_rank %>% rownames(.)
# Max 20
p1 <- plot_bc(tf_rank[1]) + plot_bc(tf_rank[2]) + plot_bc(tf_rank[3]) + plot_bc(tf_rank[4]) + plot_bc(tf_rank[5])
p2 <- plot_bc(tf_rank[6]) + plot_bc(tf_rank[7]) + plot_bc(tf_rank[8]) + plot_bc(tf_rank[9]) + plot_bc(tf_rank[10])
p3 <- plot_bc(tf_rank[11]) + plot_bc(tf_rank[12]) + plot_bc(tf_rank[13]) + plot_bc(tf_rank[14]) + plot_bc(tf_rank[15])
p4 <- plot_bc(tf_rank[16]) + plot_bc(tf_rank[17]) + plot_bc(tf_rank[18]) + plot_bc(tf_rank[19]) + plot_bc(tf_rank[20])

p5 <- plot_oc(tf_rank[1]) + plot_oc(tf_rank[2]) + plot_oc(tf_rank[3]) + plot_oc(tf_rank[4]) + plot_oc(tf_rank[5])
p6 <- plot_oc(tf_rank[6]) + plot_oc(tf_rank[7]) + plot_oc(tf_rank[8]) + plot_oc(tf_rank[9]) + plot_oc(tf_rank[10])
p7 <- plot_oc(tf_rank[11]) + plot_oc(tf_rank[12]) + plot_oc(tf_rank[13]) + plot_oc(tf_rank[14]) + plot_oc(tf_rank[15])
p8 <- plot_oc(tf_rank[16]) + plot_oc(tf_rank[17]) + plot_oc(tf_rank[18]) + plot_oc(tf_rank[19]) + plot_oc(tf_rank[20])

p9 <- plot_ic(tf_rank[1]) + plot_ic(tf_rank[2]) + plot_ic(tf_rank[3]) + plot_ic(tf_rank[4]) + plot_ic(tf_rank[5])
p10 <- plot_ic(tf_rank[6]) + plot_ic(tf_rank[7]) + plot_ic(tf_rank[8]) + plot_ic(tf_rank[9]) + plot_ic(tf_rank[10])
p11 <- plot_ic(tf_rank[11]) + plot_ic(tf_rank[12]) + plot_ic(tf_rank[13]) + plot_ic(tf_rank[14]) + plot_ic(tf_rank[15])
p12 <- plot_ic(tf_rank[16]) + plot_ic(tf_rank[17]) + plot_ic(tf_rank[18]) + plot_ic(tf_rank[19]) + plot_ic(tf_rank[20])

p1 <- grid.grabExpr(draw(p1))
p2 <- grid.grabExpr(draw(p2))
p3 <- grid.grabExpr(draw(p3))
p4 <- grid.grabExpr(draw(p4))
p5 <- grid.grabExpr(draw(p5))
p6 <- grid.grabExpr(draw(p6))
p7 <- grid.grabExpr(draw(p7))
p8 <- grid.grabExpr(draw(p8))
p9 <- grid.grabExpr(draw(p9))
p10 <- grid.grabExpr(draw(p10))
p11 <- grid.grabExpr(draw(p11))
p12 <- grid.grabExpr(draw(p12))
col_fun = colorRamp2(c(0, 0.001, 1), c('navy',"khaki1", "red"))
legend_grob1 <- grid.grabExpr(draw(Legend(col_fun = col_fun, title = "weighted\nnetwork\ncentrality")))
col_fun = colorRamp2(c(0, 0.001, 1), c('gainsboro', "grey", "black"))
legend_grob2 <- grid.grabExpr(draw(Legend(col_fun = col_fun, title = "unweighted\nnetwork\ncentrality")))

options(repr.plot.width=24, repr.plot.height=12)
plot_grid(plot_grid(textbox_grob("Main Regulatory TFs for Epidermis ", gp = gpar(fontsize = 20, col='black'), background_gp = gpar(fill = "#FFFFFF", col = NA)),
          plot_grid(plot_grid(textbox_grob("Betweenness centrality", gp = gpar(fontsize = 16, col='black'), background_gp = gpar(fill = "#FFFFFF", col = NA)),plot_grid(p1,p2,p3,p4,nrow=4),nrow=2,rel_heights = c(1, 14))
                    ,plot_grid(textbox_grob("Out-degree centrality", gp = gpar(fontsize = 16, col='black'), background_gp = gpar(fill = "#FFFFFF", col = NA)),plot_grid(p5,p6,p7,p8,nrow=4),nrow=2,rel_heights=c(1,14))
                    ,plot_grid(textbox_grob("In-degree centrality", gp = gpar(fontsize = 16, col='black'), background_gp = gpar(fill = "#FFFFFF", col = NA)),plot_grid(p9,p10,p11,p12,nrow=4),nrow=2,rel_heights=c(1,14)),nrow=1),nrow=2,rel_heights = c(1, 14)),plot_grid(legend_grob1,legend_grob2,nrow=2), rel_widths = c(19, 1))

epilrc_rank <- bc_rank[which(bc_rank$epilrc*2 > bc_rank$all),]%>% arrange(desc(epilrc))
epilrc_rank <- epilrc_rank[-c(match(rownames(atri_rank), rownames(epilrc_rank)),match(rownames(tri_rank), rownames(epilrc_rank)),match(rownames(lrc_rank), rownames(epilrc_rank))),]
epilrc_rank$GeneName <- rownames(epilrc_rank)

epilrc_rank

options(repr.plot.width=8, repr.plot.height=40)
ggplot(epilrc_rank, aes(x=reorder(GeneName, epilrc, decreasing = FALSE),  y=epilrc)) + geom_point(size=4)+ 
  labs(title="Epidermis+LRC-specific TF Prioritization",x="", y = "Combined centrality score (betweeness, out and in degree)")+
  theme_classic()+
  theme(axis.text.x = element_text(color = "black", 
                           size = 16, angle = 0),
          axis.text.y = element_text( color = "black", 
                           size = 16, angle = 0),
        plot.title = element_text(color="black", size=18, face="bold"),
        axis.title.x = element_text(color="black", size=16, face="bold"))+ coord_flip()

write.csv(epilrc_rank,"Epidermis_LRC_TF_centrality_all_transition_zscore3.csv", quote=FALSE, row.names=TRUE)

tf_rank <- epilrc_rank %>% rownames(.)
# Max 20
p1 <- plot_bc(tf_rank[1]) + plot_bc(tf_rank[2]) + plot_bc(tf_rank[3]) + plot_bc(tf_rank[4]) + plot_bc(tf_rank[5])
p2 <- plot_bc(tf_rank[6]) + plot_bc(tf_rank[7]) + plot_bc(tf_rank[8]) + plot_bc(tf_rank[9]) + plot_bc(tf_rank[10])
p3 <- plot_bc(tf_rank[11]) + plot_bc(tf_rank[12]) + plot_bc(tf_rank[13]) + plot_bc(tf_rank[14]) + plot_bc(tf_rank[15])
p4 <- plot_bc(tf_rank[16]) + plot_bc(tf_rank[17]) + plot_bc(tf_rank[18]) + plot_bc(tf_rank[19]) + plot_bc(tf_rank[20])

p5 <- plot_oc(tf_rank[1]) + plot_oc(tf_rank[2]) + plot_oc(tf_rank[3]) + plot_oc(tf_rank[4]) + plot_oc(tf_rank[5])
p6 <- plot_oc(tf_rank[6]) + plot_oc(tf_rank[7]) + plot_oc(tf_rank[8]) + plot_oc(tf_rank[9]) + plot_oc(tf_rank[10])
p7 <- plot_oc(tf_rank[11]) + plot_oc(tf_rank[12]) + plot_oc(tf_rank[13]) + plot_oc(tf_rank[14]) + plot_oc(tf_rank[15])
p8 <- plot_oc(tf_rank[16]) + plot_oc(tf_rank[17]) + plot_oc(tf_rank[18]) + plot_oc(tf_rank[19]) + plot_oc(tf_rank[20])

p9 <- plot_ic(tf_rank[1]) + plot_ic(tf_rank[2]) + plot_ic(tf_rank[3]) + plot_ic(tf_rank[4]) + plot_ic(tf_rank[5])
p10 <- plot_ic(tf_rank[6]) + plot_ic(tf_rank[7]) + plot_ic(tf_rank[8]) + plot_ic(tf_rank[9]) + plot_ic(tf_rank[10])
p11 <- plot_ic(tf_rank[11]) + plot_ic(tf_rank[12]) + plot_ic(tf_rank[13]) + plot_ic(tf_rank[14]) + plot_ic(tf_rank[15])
p12 <- plot_ic(tf_rank[16]) + plot_ic(tf_rank[17]) + plot_ic(tf_rank[18]) + plot_ic(tf_rank[19]) + plot_ic(tf_rank[20])

p1 <- grid.grabExpr(draw(p1))
p2 <- grid.grabExpr(draw(p2))
p3 <- grid.grabExpr(draw(p3))
p4 <- grid.grabExpr(draw(p4))
p5 <- grid.grabExpr(draw(p5))
p6 <- grid.grabExpr(draw(p6))
p7 <- grid.grabExpr(draw(p7))
p8 <- grid.grabExpr(draw(p8))
p9 <- grid.grabExpr(draw(p9))
p10 <- grid.grabExpr(draw(p10))
p11 <- grid.grabExpr(draw(p11))
p12 <- grid.grabExpr(draw(p12))
col_fun = colorRamp2(c(0, 0.001, 1), c('navy',"khaki1", "red"))
legend_grob1 <- grid.grabExpr(draw(Legend(col_fun = col_fun, title = "weighted\nnetwork\ncentrality")))
col_fun = colorRamp2(c(0, 0.001, 1), c('gainsboro', "grey", "black"))
legend_grob2 <- grid.grabExpr(draw(Legend(col_fun = col_fun, title = "unweighted\nnetwork\ncentrality")))

options(repr.plot.width=24, repr.plot.height=12)
plot_grid(plot_grid(textbox_grob("Main Regulatory TFs for Epidermis (includes LRC)", gp = gpar(fontsize = 20, col='black'), background_gp = gpar(fill = "#FFFFFF", col = NA)),
          plot_grid(plot_grid(textbox_grob("Betweenness centrality", gp = gpar(fontsize = 16, col='black'), background_gp = gpar(fill = "#FFFFFF", col = NA)),plot_grid(p1,p2,p3,p4,nrow=4),nrow=2,rel_heights = c(1, 14))
                    ,plot_grid(textbox_grob("Out-degree centrality", gp = gpar(fontsize = 16, col='black'), background_gp = gpar(fill = "#FFFFFF", col = NA)),plot_grid(p5,p6,p7,p8,nrow=4),nrow=2,rel_heights=c(1,14))
                    ,plot_grid(textbox_grob("In-degree centrality", gp = gpar(fontsize = 16, col='black'), background_gp = gpar(fill = "#FFFFFF", col = NA)),plot_grid(p9,p10,p11,p12,nrow=4),nrow=2,rel_heights=c(1,14)),nrow=1),nrow=2,rel_heights = c(1, 14)),plot_grid(legend_grob1,legend_grob2,nrow=2), rel_widths = c(19, 1))

stele_rank <- bc_rank[which(bc_rank$stele*2 > bc_rank$all),]%>% arrange(desc(stele))
stele_rank <- stele_rank[-c(match(rownames(per_rank), rownames(stele_rank)),match(rownames(pro_rank), rownames(stele_rank)),match(rownames(xyl_rank), rownames(stele_rank)),match(rownames(phl_rank), rownames(stele_rank))),]
stele_rank$GeneName <- rownames(stele_rank)

stele_rank

options(repr.plot.width=6, repr.plot.height=4)
ggplot(stele_rank[1:10,], aes(x=reorder(GeneName, stele, decreasing = FALSE),  y=stele)) + geom_point(size=4)+ 
  labs(title="Stele-specific TF Prioritization",x="", y = "Combined centrality score")+
  theme_classic()+
  theme(axis.text.x = element_text(color = "black", 
                           size = 16, angle = 0),
          axis.text.y = element_text( color = "black", 
                           size = 16, angle = 0),
        plot.title = element_text(color="black", size=18, face="bold"),
        axis.title.x = element_text(color="black", size=16, face="bold"))+ coord_flip()

write.csv(stele_rank,"Stele_TF_centrality_all_transition_zscore3.csv", quote=FALSE, row.names=TRUE)

tf_rank <- stele_rank %>% rownames(.)
# Max 20
p1 <- plot_bc(tf_rank[1]) + plot_bc(tf_rank[2]) + plot_bc(tf_rank[3]) + plot_bc(tf_rank[4]) + plot_bc(tf_rank[5])
p2 <- plot_bc(tf_rank[6]) + plot_bc(tf_rank[7]) + plot_bc(tf_rank[8]) + plot_bc(tf_rank[9]) + plot_bc(tf_rank[10])
p3 <- plot_bc(tf_rank[11]) + plot_bc(tf_rank[12]) + plot_bc(tf_rank[13]) + plot_bc(tf_rank[14]) + plot_bc(tf_rank[15])
p4 <- plot_bc(tf_rank[16]) + plot_bc(tf_rank[17]) + plot_bc(tf_rank[18]) + plot_bc(tf_rank[19]) + plot_bc(tf_rank[20])

p5 <- plot_oc(tf_rank[1]) + plot_oc(tf_rank[2]) + plot_oc(tf_rank[3]) + plot_oc(tf_rank[4]) + plot_oc(tf_rank[5])
p6 <- plot_oc(tf_rank[6]) + plot_oc(tf_rank[7]) + plot_oc(tf_rank[8]) + plot_oc(tf_rank[9]) + plot_oc(tf_rank[10])
p7 <- plot_oc(tf_rank[11]) + plot_oc(tf_rank[12]) + plot_oc(tf_rank[13]) + plot_oc(tf_rank[14]) + plot_oc(tf_rank[15])
p8 <- plot_oc(tf_rank[16]) + plot_oc(tf_rank[17]) + plot_oc(tf_rank[18]) + plot_oc(tf_rank[19]) + plot_oc(tf_rank[20])

p9 <- plot_ic(tf_rank[1]) + plot_ic(tf_rank[2]) + plot_ic(tf_rank[3]) + plot_ic(tf_rank[4]) + plot_ic(tf_rank[5])
p10 <- plot_ic(tf_rank[6]) + plot_ic(tf_rank[7]) + plot_ic(tf_rank[8]) + plot_ic(tf_rank[9]) + plot_ic(tf_rank[10])
p11 <- plot_ic(tf_rank[11]) + plot_ic(tf_rank[12]) + plot_ic(tf_rank[13]) + plot_ic(tf_rank[14]) + plot_ic(tf_rank[15])
p12 <- plot_ic(tf_rank[16]) + plot_ic(tf_rank[17]) + plot_ic(tf_rank[18]) + plot_ic(tf_rank[19]) + plot_ic(tf_rank[20])

p1 <- grid.grabExpr(draw(p1))
p2 <- grid.grabExpr(draw(p2))
p3 <- grid.grabExpr(draw(p3))
p4 <- grid.grabExpr(draw(p4))
p5 <- grid.grabExpr(draw(p5))
p6 <- grid.grabExpr(draw(p6))
p7 <- grid.grabExpr(draw(p7))
p8 <- grid.grabExpr(draw(p8))
p9 <- grid.grabExpr(draw(p9))
p10 <- grid.grabExpr(draw(p10))
p11 <- grid.grabExpr(draw(p11))
p12 <- grid.grabExpr(draw(p12))
col_fun = colorRamp2(c(0, 0.001, 1), c('navy',"khaki1", "red"))
legend_grob1 <- grid.grabExpr(draw(Legend(col_fun = col_fun, title = "weighted\nnetwork\ncentrality")))
col_fun = colorRamp2(c(0, 0.001, 1), c('gainsboro', "grey", "black"))
legend_grob2 <- grid.grabExpr(draw(Legend(col_fun = col_fun, title = "unweighted\nnetwork\ncentrality")))

options(repr.plot.width=24, repr.plot.height=12)
plot_grid(plot_grid(textbox_grob("Main Regulatory TFs for Stele", gp = gpar(fontsize = 20, col='black'), background_gp = gpar(fill = "#FFFFFF", col = NA)),
          plot_grid(plot_grid(textbox_grob("Betweenness centrality", gp = gpar(fontsize = 16, col='black'), background_gp = gpar(fill = "#FFFFFF", col = NA)),plot_grid(p1,p2,p3,p4,nrow=4),nrow=2,rel_heights = c(1, 14))
                    ,plot_grid(textbox_grob("Out-degree centrality", gp = gpar(fontsize = 16, col='black'), background_gp = gpar(fill = "#FFFFFF", col = NA)),plot_grid(p5,p6,p7,p8,nrow=4),nrow=2,rel_heights=c(1,14))
                    ,plot_grid(textbox_grob("In-degree centrality", gp = gpar(fontsize = 16, col='black'), background_gp = gpar(fill = "#FFFFFF", col = NA)),plot_grid(p9,p10,p11,p12,nrow=4),nrow=2,rel_heights=c(1,14)),nrow=1),nrow=2,rel_heights = c(1, 14)),plot_grid(legend_grob1,legend_grob2,nrow=2), rel_widths = c(19, 1))

rc_rank <- bc_rank[which(bc_rank$rc*2 > bc_rank$all),]%>% arrange(desc(rc))
rc_rank <- rc_rank[-c(match(rownames(lrc_rank), rownames(rc_rank)),match(rownames(col_rank), rownames(rc_rank))),]
rc_rank$GeneName <- rownames(rc_rank)

rc_rank

rc_rank[1:50,]

options(repr.plot.width=6, repr.plot.height=4)
ggplot(rc_rank[1:10,], aes(x=reorder(GeneName, rc, decreasing = FALSE),  y=rc)) + geom_point(size=4)+ 
  labs(title="Root Cap-specific TF Prioritization",x="", y = "Combined centrality score")+
  theme_classic()+
  theme(axis.text.x = element_text(color = "black", 
                           size = 16, angle = 0),
          axis.text.y = element_text( color = "black", 
                           size = 16, angle = 0),
        plot.title = element_text(color="black", size=18, face="bold"),
        axis.title.x = element_text(color="black", size=16, face="bold"))+ coord_flip()

write.csv(rc_rank,"Root_Cap_TF_centrality_all_transition_zscore3.csv", quote=FALSE, row.names=TRUE)

tf_rank <- rc_rank %>% rownames(.)
# Max 20
p1 <- plot_bc(tf_rank[1]) + plot_bc(tf_rank[2]) + plot_bc(tf_rank[3]) + plot_bc(tf_rank[4]) + plot_bc(tf_rank[5])
p2 <- plot_bc(tf_rank[6]) + plot_bc(tf_rank[7]) + plot_bc(tf_rank[8]) + plot_bc(tf_rank[9]) + plot_bc(tf_rank[10])
p3 <- plot_bc(tf_rank[11]) + plot_bc(tf_rank[12]) + plot_bc(tf_rank[13]) + plot_bc(tf_rank[14]) + plot_bc(tf_rank[15])
p4 <- plot_bc(tf_rank[16]) + plot_bc(tf_rank[17]) + plot_bc(tf_rank[18]) + plot_bc(tf_rank[19]) + plot_bc(tf_rank[20])

p5 <- plot_oc(tf_rank[1]) + plot_oc(tf_rank[2]) + plot_oc(tf_rank[3]) + plot_oc(tf_rank[4]) + plot_oc(tf_rank[5])
p6 <- plot_oc(tf_rank[6]) + plot_oc(tf_rank[7]) + plot_oc(tf_rank[8]) + plot_oc(tf_rank[9]) + plot_oc(tf_rank[10])
p7 <- plot_oc(tf_rank[11]) + plot_oc(tf_rank[12]) + plot_oc(tf_rank[13]) + plot_oc(tf_rank[14]) + plot_oc(tf_rank[15])
p8 <- plot_oc(tf_rank[16]) + plot_oc(tf_rank[17]) + plot_oc(tf_rank[18]) + plot_oc(tf_rank[19]) + plot_oc(tf_rank[20])

p9 <- plot_ic(tf_rank[1]) + plot_ic(tf_rank[2]) + plot_ic(tf_rank[3]) + plot_ic(tf_rank[4]) + plot_ic(tf_rank[5])
p10 <- plot_ic(tf_rank[6]) + plot_ic(tf_rank[7]) + plot_ic(tf_rank[8]) + plot_ic(tf_rank[9]) + plot_ic(tf_rank[10])
p11 <- plot_ic(tf_rank[11]) + plot_ic(tf_rank[12]) + plot_ic(tf_rank[13]) + plot_ic(tf_rank[14]) + plot_ic(tf_rank[15])
p12 <- plot_ic(tf_rank[16]) + plot_ic(tf_rank[17]) + plot_ic(tf_rank[18]) + plot_ic(tf_rank[19]) + plot_ic(tf_rank[20])

p1 <- grid.grabExpr(draw(p1))
p2 <- grid.grabExpr(draw(p2))
p3 <- grid.grabExpr(draw(p3))
p4 <- grid.grabExpr(draw(p4))
p5 <- grid.grabExpr(draw(p5))
p6 <- grid.grabExpr(draw(p6))
p7 <- grid.grabExpr(draw(p7))
p8 <- grid.grabExpr(draw(p8))
p9 <- grid.grabExpr(draw(p9))
p10 <- grid.grabExpr(draw(p10))
p11 <- grid.grabExpr(draw(p11))
p12 <- grid.grabExpr(draw(p12))
col_fun = colorRamp2(c(0, 0.001, 1), c('navy',"khaki1", "red"))
legend_grob1 <- grid.grabExpr(draw(Legend(col_fun = col_fun, title = "weighted\nnetwork\ncentrality")))
col_fun = colorRamp2(c(0, 0.001, 1), c('gainsboro', "grey", "black"))
legend_grob2 <- grid.grabExpr(draw(Legend(col_fun = col_fun, title = "unweighted\nnetwork\ncentrality")))

options(repr.plot.width=24, repr.plot.height=12)
plot_grid(plot_grid(textbox_grob("Main Regulatory TFs for Root Cap", gp = gpar(fontsize = 20, col='black'), background_gp = gpar(fill = "#FFFFFF", col = NA)),
          plot_grid(plot_grid(textbox_grob("Betweenness centrality", gp = gpar(fontsize = 16, col='black'), background_gp = gpar(fill = "#FFFFFF", col = NA)),plot_grid(p1,p2,p3,p4,nrow=4),nrow=2,rel_heights = c(1, 14))
                    ,plot_grid(textbox_grob("Out-degree centrality", gp = gpar(fontsize = 16, col='black'), background_gp = gpar(fill = "#FFFFFF", col = NA)),plot_grid(p5,p6,p7,p8,nrow=4),nrow=2,rel_heights=c(1,14))
                    ,plot_grid(textbox_grob("In-degree centrality", gp = gpar(fontsize = 16, col='black'), background_gp = gpar(fill = "#FFFFFF", col = NA)),plot_grid(p9,p10,p11,p12,nrow=4),nrow=2,rel_heights=c(1,14)),nrow=1),nrow=2,rel_heights = c(1, 14)),plot_grid(legend_grob1,legend_grob2,nrow=2), rel_widths = c(19, 1))

	used	(Mb)	gc trigger	(Mb)	max used	(Mb)
Ncells	625293	33.4	1361470	72.8	1361470	72.8
Vcells	1159120	8.9	8388608	64.0	1802279	13.8

	X	tf_occurance	atri_degree_centrality	atri_out_centrality	atri_in_centrality	atri_betweenness_centrality	atri_closeness_centrality	atri_eigenvector_centrality	tri_degree_centrality	tri_out_centrality	...	phl_in_centrality	phl_betweenness_centrality	phl_closeness_centrality	phl_eigenvector_centrality	col_degree_centrality	col_out_centrality	col_in_centrality	col_betweenness_centrality	col_closeness_centrality	col_eigenvector_centrality
	<chr>	<dbl>	<dbl>	<dbl>	<dbl>	<dbl>	<dbl>	<dbl>	<dbl>	<dbl>	...	<dbl>	<dbl>	<dbl>	<dbl>	<dbl>	<dbl>	<dbl>	<dbl>	<dbl>	<dbl>
1	H2AXA	10	0.4023372	0.1969950	0.2053422	0.04084567	0.002893454	0.06930489	0.3385580	0.10658307	...	0.20268757	0.1068750973	0.0011064411	0.04917093	0.3140625	0.2375000	0.0765625	0.4845462	0.003069275	0.05594356
2	HB-2	10	1.0217028	0.2904841	0.7312187	0.88444230	0.003620554	0.10318643	0.2915361	0.03134796	...	0.03471445	0.0006063604	0.0011005165	0.05188664	1.1687500	0.6250000	0.5437500	0.9751589	0.004531334	0.10857313
3	GATA19	10	1.0200334	0.6327212	0.3873122	0.98922117	0.003523433	0.10154627	0.2539185	0.10344828	...	0.13101904	0.4793435741	0.0009728508	0.04139889	0.2281250	0.0609375	0.1671875	0.9294601	0.003688327	0.04870697
4	HTA10	10	0.4590985	0.2237062	0.2353923	0.56394716	0.002961137	0.07343545	0.3855799	0.12852665	...	0.28891377	0.4712185961	0.0011567435	0.05680820	0.2781250	0.1593750	0.1187500	0.6063503	0.002957289	0.05189729
5	GAMMA-H2AX	10	0.5876461	0.3756260	0.2120200	0.87042786	0.003067845	0.07959582	0.8275862	0.50156740	...	0.27547592	0.0141747724	0.0012103003	0.05417356	0.3062500	0.1546875	0.1515625	0.9231123	0.003441618	0.05485739
6	GRP2B	9	0.3405676	0.1903172	0.1502504	0.90321104	0.002998468	0.06278893	0.0000000	0.00000000	...	0.20716685	0.9618758756	0.0011917313	0.08047918	0.5062500	0.4484375	0.0578125	0.4854582	0.003975758	0.08353683

	X	tf_occurance	atri_degree_centrality	atri_out_centrality	atri_in_centrality	atri_betweenness_centrality	atri_closeness_centrality	atri_eigenvector_centrality	tri_degree_centrality	tri_out_centrality	...	phl_in_centrality	phl_betweenness_centrality	phl_closeness_centrality	phl_eigenvector_centrality	col_degree_centrality	col_out_centrality	col_in_centrality	col_betweenness_centrality	col_closeness_centrality	col_eigenvector_centrality
	<chr>	<dbl>	<dbl>	<dbl>	<dbl>	<dbl>	<dbl>	<dbl>	<dbl>	<dbl>	...	<dbl>	<dbl>	<dbl>	<dbl>	<dbl>	<dbl>	<dbl>	<dbl>	<dbl>	<dbl>
1	H2AXA	7	0.1624473	0.08227848	0.08016878	0.01555294	0.001103735	0.05928981	0.14465409	0.058700210	...	0.21161049	0.000000000	0.0010778420	0.05035964	0.35972222	0.086111111	0.27361111	0.07714225	0.0007635986	0.08372615
2	BME3	7	0.2616034	0.22995781	0.03164557	0.03877753	0.001300033	0.07647016	0.24528302	0.140461216	...	0.03183521	0.641672112	0.0006751422	0.02687088	0.04722222	0.005555556	0.04166667	0.00000000	0.0006292341	0.02312157
3	HAT1	7	0.3080169	0.15611814	0.15189873	0.83841803	0.001372900	0.08713217	0.07337526	0.008385744	...	0.12921348	0.856219828	0.0007598532	0.04521189	0.13194444	0.034722222	0.09722222	0.80140434	0.0006157478	0.04761614
4	MYB12	6	0.4324895	0.40717300	0.02531646	0.18455232	0.001544689	0.09978931	0.03563941	0.016771488	...	0.00000000	0.000000000	0.0000000000	0.00000000	0.00000000	0.000000000	0.00000000	0.00000000	0.0000000000	0.00000000
5	GRP2B	6	0.2531646	0.14556962	0.10759494	0.04942418	0.001365611	0.07830556	0.00000000	0.000000000	...	0.05430712	0.031413594	0.0010275356	0.08967945	0.08333333	0.029166667	0.05416667	0.42131819	0.0004438084	0.03616124
6	GATA15	6	0.3797468	0.19620253	0.18354430	0.06227866	0.001610540	0.09551755	0.23270440	0.041928721	...	0.02996255	0.002817772	0.0005445657	0.01997458	0.00000000	0.000000000	0.00000000	0.00000000	0.0000000000	0.00000000

	X	tf_occurance	atri_degree_centrality	atri_out_centrality	atri_in_centrality	atri_betweenness_centrality	atri_closeness_centrality	atri_eigenvector_centrality	tri_degree_centrality	tri_out_centrality	...	phl_in_centrality	phl_betweenness_centrality	phl_closeness_centrality	phl_eigenvector_centrality	col_degree_centrality	col_out_centrality	col_in_centrality	col_betweenness_centrality	col_closeness_centrality	col_eigenvector_centrality
	<chr>	<dbl>	<dbl>	<dbl>	<dbl>	<dbl>	<dbl>	<dbl>	<dbl>	<dbl>	...	<dbl>	<dbl>	<dbl>	<dbl>	<dbl>	<dbl>	<dbl>	<dbl>	<dbl>	<dbl>
1	AT4G30410	6	0.7468354	0.60506329	0.1417722	0.856499390	0.0004675479	0.15472441	0.8629550	0.53319058	...	0.00000000	0.0000000	0.0000000000	0.00000000	0.00000000	0.000000000	0.00000000	0.0000000	0.0000000000	0.00000000
2	ATCTH	6	0.3088608	0.05569620	0.2531646	0.002081861	0.0004639463	0.10061684	0.1605996	0.02569593	...	0.00000000	0.0000000	0.0000000000	0.00000000	0.35532995	0.299492386	0.05583756	0.7610230	0.0010902563	0.08569533
3	IAA14	5	0.3392405	0.25569620	0.0835443	0.009760329	0.0003937151	0.10527274	0.2312634	0.13918630	...	0.00000000	0.0000000	0.0000000000	0.00000000	0.01649746	0.001269036	0.01522843	0.0000000	0.0007651091	0.01138349
4	BIM1	5	0.5316456	0.31645570	0.2151899	0.011379554	0.0004630263	0.13175744	0.1070664	0.06209850	...	0.07387141	0.0000000	0.0003685896	0.04403974	0.27918782	0.189086294	0.09010152	0.9445656	0.0012919391	0.07964298
5	GATA4	5	0.2911392	0.04556962	0.2455696	0.001702757	0.0005350022	0.09326911	0.6295503	0.02783726	...	0.00000000	0.0000000	0.0000000000	0.00000000	0.00000000	0.000000000	0.00000000	0.0000000	0.0000000000	0.00000000
6	HAT2	5	0.0000000	0.00000000	0.0000000	0.000000000	0.0000000000	0.00000000	0.0000000	0.00000000	...	0.07250342	0.9518112	0.0003552297	0.05053935	0.00000000	0.000000000	0.00000000	0.0000000	0.0000000000	0.00000000

	X	tf_occurance	atri_degree_centrality	atri_out_centrality	atri_in_centrality	atri_betweenness_centrality	atri_closeness_centrality	atri_eigenvector_centrality	tri_degree_centrality	tri_out_centrality	...	phl_in_centrality	phl_betweenness_centrality	phl_closeness_centrality	phl_eigenvector_centrality	col_degree_centrality	col_out_centrality	col_in_centrality	col_betweenness_centrality	col_closeness_centrality	col_eigenvector_centrality
	<chr>	<dbl>	<dbl>	<dbl>	<dbl>	<dbl>	<dbl>	<dbl>	<dbl>	<dbl>	...	<dbl>	<dbl>	<dbl>	<dbl>	<dbl>	<dbl>	<dbl>	<dbl>	<dbl>	<dbl>
1	RD26	7	0.6697248	0.43425076	0.23547401	0.145531979	0.0006467123	0.13050671	0.16901408	0.118309859	...	0.004237288	0.004712759	0.0003913220	0.04770843	0.01111111	0.003174603	0.007936508	0.000000000	0.000717119	0.006922739
2	HAT22	6	0.1345566	0.09174312	0.04281346	0.000000000	0.0004387974	0.04661033	0.53802817	0.400000000	...	0.032485876	0.000000000	0.0002677743	0.02895724	0.08412698	0.019047619	0.065079365	0.001249148	0.001133165	0.037746444
3	MYB7	6	0.8379205	0.45259939	0.38532110	0.947468153	0.0007034644	0.13880090	0.00000000	0.000000000	...	0.000000000	0.000000000	0.0000000000	0.00000000	0.00000000	0.000000000	0.000000000	0.000000000	0.000000000	0.000000000
4	KNAT7	6	0.3822630	0.33639144	0.04587156	0.289581809	0.0006074321	0.09458828	0.01126761	0.005633803	...	0.000000000	0.000000000	0.0000000000	0.00000000	0.71904762	0.284126984	0.434920635	0.852792288	0.001607269	0.122205425
5	ATHB13	6	0.3363914	0.03058104	0.30581040	0.018348624	0.0007202355	0.08535107	0.06760563	0.011267606	...	0.000000000	0.000000000	0.0000000000	0.00000000	0.01904762	0.003174603	0.015873016	0.056128397	0.001144112	0.010720118
6	RAV1	6	0.4006116	0.22935780	0.17125382	0.004493349	0.0005668635	0.10055625	0.10704225	0.084507042	...	0.056497175	0.952626679	0.0002896641	0.05385743	0.14285714	0.092063492	0.050793651	0.450818886	0.001049787	0.047376558

Load TF list¶

Load GRN centrality scores for each transition (LIMA) (GeneSys)¶

Load GRN centrality scores for each transition (Linear Interaction Matrix, LIMA) (scRNAseq)¶

GO terms (Prepare GS)¶

R50¶

LIMA (GeneSys)¶

LIMA (scRNAseq)¶

CellOracle¶

Differential Expression Analysis¶

Plot r50 with ranked by different centrality¶

R50 for Tissue specific TFs¶

Stele (xylem, phloem, procambium)¶

Epidermis¶

R50 for celltype specific¶

Xylem¶

Trichoblast¶

Plotting¶

Rank by betweenness centrality, out-degree centrality, in-degree centrality and degree centrality¶

Atrichoblast¶

Trichoblast¶

Cortex¶

Endodermis¶

Pericycle¶

Procambium¶

Phloem¶

Xylem¶

LRC¶

Columella¶

Add tissue level¶

Ground Tissue¶

Epidermis¶

Epidermis + LRC¶

Stele¶

Root Cap¶

Methods	R50
<chr>	<dbl>
LIMA degree centrality (GeneSys)	314.000
LIMA out centrality (GeneSys)	327.000
LIMA in centrality (GeneSys)	375.000
LIMA betweenness centrality (GeneSys)	335.000
LIMA eigenvector centrality (GeneSys)	432.000
LIMA degree centrality (scRNAseq)	365.000
LIMA out centrality (scRNAseq)	330.000
LIMA in centrality (scRNAseq)	444.000
LIMA betweenness centrality (scRNAseq)	447.000
LIMA eigenvector centrality (scRNAseq)	536.000
CellOracle degree centrality (scRNAseq)	470.000
CellOracle out centrality (scRNAseq)	513.000
CellOracle in centrality (scRNAseq)	517.000
CellOracle betweenness centrality (scRNAseq)	410.000
CellOracle eigenvector centrality (scRNAseq)	453.000
CellOracle degree centrality (GeneSys)	434.000
CellOracle out centrality (GeneSys)	466.000
CellOracle in centrality (GeneSys)	312.000
CellOracle betweenness centrality (GeneSys)	335.000
CellOracle eigenvector centrality (GeneSys)	445.000
DE myAUC rank (scRNAseq)	380.000
DE avg diff rank (scRNAseq)	403.000
DE myAUC rank (GeneSys)	310.000
DE avg diff rank (GeneSys)	331.000
Permutation	743.406

A data.frame: 6 x 62
	X	tf_occurance	atri_degree_centrality	atri_out_centrality	atri_in_centrality	atri_betweenness_centrality	atri_closeness_centrality	atri_eigenvector_centrality	tri_degree_centrality	tri_out_centrality	...	phl_in_centrality	phl_betweenness_centrality	phl_closeness_centrality	phl_eigenvector_centrality	col_degree_centrality	col_out_centrality	col_in_centrality	col_betweenness_centrality	col_closeness_centrality	col_eigenvector_centrality
	<chr>	<dbl>	<dbl>	<dbl>	<dbl>	<dbl>	<dbl>	<dbl>	<dbl>	<dbl>	...	<dbl>	<dbl>	<dbl>	<dbl>	<dbl>	<dbl>	<dbl>	<dbl>	<dbl>	<dbl>
1	HB-1	8	0.2195122	0.1324042	0.08710801	0.0002802076	0.0004585878	0.06551537	0.10238908	0.027303754	...	0.02790015	0.000000e+00	0.0003329542	0.04644030	0.41637631	0.109756098	0.30662021	0.3944822470	0.0015520636	0.09681453
2	WLIM1	7	0.5296167	0.2891986	0.24041812	0.9633171706	0.0005781904	0.10994336	0.10238908	0.088737201	...	0.03817915	1.390688e-03	0.0004631928	0.05127801	0.00000000	0.000000000	0.00000000	0.0000000000	0.0000000000	0.00000000
3	RD26	7	0.7386760	0.5017422	0.23693380	0.9363441437	0.0006638061	0.13094847	0.27303754	0.187713311	...	0.04405286	1.943509e-05	0.0003849459	0.04142436	0.01916376	0.006968641	0.01219512	0.0000000000	0.0009426614	0.01042103
4	MYB7	7	0.6236934	0.4250871	0.19860627	0.9121244609	0.0006516830	0.11882806	0.02730375	0.003412969	...	0.00000000	0.000000e+00	0.0000000000	0.00000000	0.00000000	0.000000000	0.00000000	0.0000000000	0.0000000000	0.00000000
5	KNAT7	7	0.5783972	0.3972125	0.18118467	0.2351063571	0.0006330818	0.11903518	0.00000000	0.000000000	...	0.00000000	0.000000e+00	0.0000000000	0.00000000	0.62020906	0.268292683	0.35191638	0.8681339730	0.0013507501	0.11408248
6	BLH3	6	0.1428571	0.1149826	0.02787456	0.0151799420	0.0004492747	0.04630695	0.01365188	0.010238908	...	0.00000000	0.000000e+00	0.0000000000	0.00000000	0.11672474	0.085365854	0.03135889	0.0004317395	0.0012679092	0.04344304

A data.frame: 6 x 302
	GeneID	TF	atri_degree_centrality_1	atri_out_centrality_1	atri_in_centrality_1	atri_betweenness_centrality_1	atri_closeness_centrality_1	atri_eigenvector_centrality_1	tri_degree_centrality_1	tri_out_centrality_1	...	phl_in_centrality_5	phl_betweenness_centrality_5	phl_closeness_centrality_5	phl_eigenvector_centrality_5	col_degree_centrality_5	col_out_centrality_5	col_in_centrality_5	col_betweenness_centrality_5	col_closeness_centrality_5	col_eigenvector_centrality_5
	<chr>	<chr>	<chr>	<chr>	<chr>	<chr>	<chr>	<chr>	<chr>	<chr>	...	<chr>	<chr>	<chr>	<chr>	<chr>	<chr>	<chr>	<chr>	<chr>	<chr>
1	AT1G08880	H2AXA	0.280558789289872	0.276995305164319	0.237911025145068	0.040883332588188	0.77982289908458	0.638016012368832	0.286472148541114	0.147186147186147	...	0.0161527165932452	0	0.5331869128049	0.0917789928495729	0.0385462555066079	0.0127931769722815	0.061965811965812	0.00594437414013402	0.644609523078635	0.182176552700344
2	AT4G16780	HB-2	0.712456344586729	0.408450704225352	0.847195357833656	0.88525785755801	0.975785470219646	0.949927136438232	0.246684350132626	0.0432900432900433	...	0	0	0	0	0	0	0	0	0	0
3	AT4G36620	GATA19	0.711292200232829	0.889671361502347	0.448742746615087	0.990133343765369	0.949610304448685	0.934827955774441	0.214854111405836	0.142857142857143	...	0	0	0	0	0	0	0	0	0	0
4	AT1G51060	HTA10	0.320139697322468	0.314553990610329	0.272727272727273	0.564467183797559	0.798064297939221	0.676041679289166	0.326259946949602	0.177489177489178	...	0.0220264317180617	0	0.568027911267958	0.120142431316983	0.00770925110132159	0.00213219616204691	0.0128205128205128	0	0.483551661993331	0.0541577050018262
5	AT1G54690	GAMMA-H2AX	0.409778812572759	0.528169014084507	0.245647969052224	0.871230495819734	0.826823493068053	0.732753574021271	0.70026525198939	0.692640692640693	...	0	0	0	0	0	0	0	0	0	0
6	AT2G21060	GRP2B	0.237485448195576	0.267605633802817	0.174081237911025	0.904043903965664	0.808125480825825	0.578030536398732	0	0	...	0.110132158590308	0.00384910102502027	0.807871207844227	0.612066867831852	0.00991189427312775	0.00426439232409382	0.014957264957265	0	0.459717836582851	0.0640027480741151

A data.frame: 6 x 305
	GeneID	TF	atri_degree_centrality_1	atri_out_centrality_1	atri_in_centrality_1	atri_betweenness_centrality_1	atri_closeness_centrality_1	atri_eigenvector_centrality_1	tri_degree_centrality_1	tri_out_centrality_1	...	phl_eigenvector_centrality_5	col_degree_centrality_5	col_out_centrality_5	col_in_centrality_5	col_betweenness_centrality_5	col_closeness_centrality_5	col_eigenvector_centrality_5	combined_score	celltype_specificity	weighted_score
	<chr>	<chr>	<chr>	<chr>	<chr>	<chr>	<chr>	<chr>	<chr>	<chr>	...	<chr>	<chr>	<chr>	<chr>	<chr>	<chr>	<chr>	<dbl>	<dbl>	<dbl>
1	AT5G24800	BZIP9	0	0	0	0	0	0	0	0	...	0.965592719930867	0	0	0	0	0	0	0.7720647	0.5822222	1.354287
2	AT3G43430	AT3G43430	0	0	0	0	0	0	0	0	...	0.390917063818791	0	0	0	0	0	0	0.9221986	0.4222222	1.344421
3	AT2G45050	GATA2	0.828870779976717	0.995305164319249	0.55705996131528	0.578715294854026	0.984309713319463	0.979686642330802	0.838196286472148	0.961038961038961	...	0	0	0	0	0	0	0	0.8763536	0.4311111	1.307465
4	AT5G13910	LEP	0	0	0	0	0	0	0.10079575596817	0.155844155844156	...	0	0	0	0	0	0	0	1.0000000	0.3066667	1.306667
5	AT1G66230	MYB20	0	0	0	0	0	0	0	0	...	0.0606489225228414	0	0	0	0	0	0	0.7950683	0.4844444	1.279513
6	AT1G07640	OBP2	0	0	0	0	0	0	0	0	...	0.982494065594694	0	0	0	0	0	0	0.5664295	0.6800000	1.246430

A data.frame: 6 x 62
	X	tf_occurance	atri_degree_centrality	atri_out_centrality	atri_in_centrality	atri_betweenness_centrality	atri_closeness_centrality	atri_eigenvector_centrality	tri_degree_centrality	tri_out_centrality	...	phl_in_centrality	phl_betweenness_centrality	phl_closeness_centrality	phl_eigenvector_centrality	col_degree_centrality	col_out_centrality	col_in_centrality	col_betweenness_centrality	col_closeness_centrality	col_eigenvector_centrality
	<chr>	<dbl>	<dbl>	<dbl>	<dbl>	<dbl>	<dbl>	<dbl>	<dbl>	<dbl>	...	<dbl>	<dbl>	<dbl>	<dbl>	<dbl>	<dbl>	<dbl>	<dbl>	<dbl>	<dbl>
1	BBX29	10	0.2473118	0.05376344	0.19354839	0.0024418167	0.0002779949	0.07634188	0.11006289	0.006289308	...	0.08914100	0.0360194910	1.067095e-04	0.12047289	0.4209184	0.1403061	0.28061224	4.642662e-04	0.0002109170	0.06890260
2	ATWHY2	10	0.3252688	0.29032258	0.03494624	0.0008767354	0.0003126635	0.08784270	0.10377358	0.081761006	...	0.03079417	0.0002762634	8.076810e-05	0.04041359	0.1275510	0.0752551	0.05229592	3.810241e-03	0.0001721267	0.03419012
3	HAT1	10	0.2795699	0.11290323	0.16666667	0.0092383271	0.0002821786	0.07954667	0.07861635	0.003144654	...	0.12965964	0.0188043318	1.094541e-04	0.13846171	0.3558673	0.1160714	0.23979592	2.687857e-04	0.0002097424	0.06503697
4	HMGB5	10	0.7204301	0.52150538	0.19892473	0.0658058720	0.0002974142	0.13050351	0.34591195	0.248427673	...	0.11669368	0.0003288851	1.050922e-04	0.06751047	0.4311224	0.1568878	0.27423469	3.502359e-04	0.0002140065	0.07247327
5	AT4G25210	10	0.1344086	0.03763441	0.09677419	0.0001956352	0.0002611916	0.05197195	0.04402516	0.009433962	...	0.04862237	0.0003183607	8.684771e-05	0.05319400	0.1441327	0.0255102	0.11862245	6.516016e-06	0.0001854091	0.03889094
6	HB-2	10	0.5645161	0.09408602	0.47043011	0.0330406052	0.0003468205	0.11220601	0.37735849	0.015723270	...	0.35818476	0.2002541624	1.248236e-04	0.15678872	1.1607143	0.4107143	0.75000000	7.682383e-03	0.0002419695	0.10843956